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ABSTRACT 


Underintegrated methods are investigated with respect to their 
stability and convergence properties. The focus was on 
identifying regions where they work and regions where techniques 
such as hourglass viscosity and hourglass control can be used. 
Results obtained show that underintegrated methods typically lead 
to finite element stiffness with spurious modes in the solution. 
However, problems exist (scalar elliptic boundary value problems) 
where underintegrated with hourglass control yield convergent 
solutions. Also, stress averaging in underintegrated stiffness 
calculations does not necessarily lead to stable or convergent 
stress states. 





0 . INTRODUCTION 


One of the most important and widespread numerical procedures used 
in contemporary finite element analysis of nonlinear problems in structural 
mechanics is the use of so-called reduced or underintegration. Promoted 
strongly in the late seventies and early eighties as a means for dramatically 
reducing computational times in large-scale calculations, the use of under- 
integrated finite element methods has become common practice in a large 
majority of all nonlinear calculations. 

A question of overriding importance that has perplexed many users 
of underintegrated finite element techniques for some years is whether 
or not these underintegrated methods are really satisfactory. It is 
known that underintegrated methods are frequently unstable, but these 
instabilities can be dampened out by the use of various types of "hourglass 

viscosity" or "hourglass control." It is also known that in many cases 

the underintegrated solutions can seem to converge at a rate equal to 

that of the fully integrated solutions. What are the true properties 

of underintegrated methods? When do they work? What criteria must 

hold in order that they can be used with confidence? 

These are the, questions that were addressed in the research project 
reported in this document. This final report summarizes the results 

of a two-year research project, supported by the NASA Lewis Research 
Center and carried out by Professor J. Tinsley Oden and his students 
at The University of Texas. Some of the principal conclusions of the 
work are listed as follows: 

1) Underintegration of finite element stiffnesses generally leads 
to the introduction of spurious modes in the finite element solution. 
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These spurious modes arise from two distinctly different mechanisms: under- 
integration of constraint terms, which gives rise to checkerboard instabili- 
ties and the underintegration of primary stiffness terms, which gives 
rise to hourglass instabilities. 

2) The spurious modes actually arise from expanded kernels of constraint 
operator and the governing differential operator. For example, an improperly 
underintegrated stiffness matrix will be ranked deficient 'and these ranked 
deficiencies correspond to additional modes supplied to the rigid body 
modes that appropriately belong to the kernel of this operator. In 
a similar fashion, underintegration of constraint terms leads to checkerboard 
modes, which belong to an expanded kernel of the constraint operator. 

3) There is a significant class of problems in which, with appropriate 
filtering, can be shown that an underintegrated solution with hourglass 
control can yield very satisfactory answers, and produce a finite element 
method which has the same rate of convergence as the fully integrated 
method. The fact that this does indeed hold has been rigorously proved 
in the enclosed document for a class of scalar elliptic boundary value 
problems . 

4) Unfortunately, underintegrated with hourglass control does not 
work uniformly on all linear or nonlinear problems, and it can lead 
to solutions which, while looking reasonable to the unsuspecting eye, 
may be grossly in error. The success of underintegrated methods seems 
to depend strongly on the regularity of the solution. Underintegration 
seems to work well in the presence of smooth solutions. 

5) Host of the better known and often used underintegrated methods 
for constrained problems are actually unstable, but the instabilities 
are subtle and may be manifested only in cases in which irregular meshes 



are used or in which there are irregularities in the data. In general, 
these unstable methods should be avoided in code development. 

6) For the underintegration of constraints, such as those occurring 
involved with the incompressiblity condition and Stokes problems are incom- 
pressible elasticity or incompressible plasticity, a necessary condition 
for the numerical stability of underintegrated methods is the satisfaction 
of a specific LBB condition. Some excellent underintegf ated elements 
which satisfy this condition are discussed in the report. Stress averaging 
in underintegrated stiffness calculations does not necessarily lead to 
a stable or convergent stress. 

0.1. Major Publications and Presentations 


A number of significant papers and reports were published during 
the contract period. These are listed as follows: 

Oden, J.T., "Penalty Method and Reduced Integration for the Analysis 
of Fluids," Proceedings, Symposium on Penalty Finite Element Methods in 
Mechanics , ASME Winter Annual Meeting, November 14-19, 1982, Phoenix, 
AZ. * . . 


Oden, J.T. and O.-P. Jacquotte, "Stability of Some Mixed Finite 
Element Methods for Stokesian Flows," Computer Methods in Applied Mechanics 
and Engineering , 1984, Vol. 43, No. 2, pp. 231-248. 

Kikuchi, N. , Oden, J.T., and Song, Y.J. , "Convergence of Modified 
Penalty Methods and Smoothing Schemes of Pressure for Stokes Flow Problems," 
Finite Elements in Fluid Dynamics, Vol. V , John Wiley & Sons, Ltd., 
London, 1984. 

Oden, J.T. and Jacquotte, O.-P., "Stable and Unstable RIP/Perturbed 
Lagrangian Methods for Two-Dimensional Viscous Flow Problems," Finite Elements 
in Fluid Dynamics, Vol. V, John Wiley & Sons, Ltd., London, 1984, pp. 
127-146. 

Endo, T., Oden, J.T. , Becker, E. and T. Miller, "A Numerical Analysis 
of Contact and Limit-Point Behavior in a Class of Problems of Finite 
Elastic Deformation," Computers and Structures, 1984, Vol. 18, No. 5, 
pp. 899-910. 

Jacquotte, O.-P. and Oden, J.T. Analysis of Hourglass Instabilities 
and Control in Underintegrated Finite Element Methods," Computer Methods 
in Applied Mechanics and Engineering , 1984, Vol. 44, pp. 339-363. 

Jacquotte, O.-P. and Oden, J.T. "Analysis and Treatment of Hourglass 
Instabilities in Underintegrated Finite Element Methods," Proceedings, Sympo- 
sium on Innovative Methods for Nonlinear Mechanics , ASME Winter Annual 
Meeting, December 1215, 1984, New Orleans, LA. 
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Jacquotte, O.-P. , "Stability , Accuracy, and Efficiency of Some Underin- 
tegrated Methods in Finite Element Computations," Computer Methods in 
Applied Mechanics and Engineering , (to appear). 

Oden, J.T., Jacquotte, O.-P. and Becker, E.B. ,"Numerical Control of 
Hourglass Instability," Computers and Structures , (to appear). 

Oden, J.T. and Jacquotte, O.-P. j"Convergence and Stability of Underin- 
tegrated Finite Element Methods," To appear in Proceedings , ASCE/ASME 
Mechanics Meeting , June 24-26, 1985 at Albuquerque, NM. 

0.2 Dissertations : 


Jacquotte, Olivier-P. , "Underintegration in Finite Element Methods," 
Ph.D. thesis. University of Texas, Austin, Texas, 1985. 

0.3 Oral Presentations 


There were three oral presentations during the research period. They 
are as follows: 

Oden, J.T., "Stability and Convergence of Underintegrated Finite Element 
Approximations," Presented at the NASA-LeRC/INDUSTRY/UNIVERSITY Workshop 
on Nonlinear Analyses for Engine Structures, April 19-20, 1983, in Cleveland, 

OH. . ■ * . . . ... 


Jacquotte, O.-P., "Analysis and Treatment of Hourglass Instabilities 
in Underintegrated Finite Element Methods," Presented at the ASME Winter 
Annual Meeting December 12-15, 1984 in New Orleans, LA. 

Oden, J.T., "Convergence and Stability of Underintegrated Finite Element 
Methods," Presented at the ASCE/ASME Mechanics Meeting June 24-26, 1985 
in Albuquerque, NM. 

0.4 Personnel . 

The following individuals worked on technical aspects of the project 
during the report period: 

Prof. J.T. Oden, Principal Investigator 

j 

Mr. O.-P. Jacquotte, Graduate Research Assistant 

Mssrs. Lin, Martins, Strouboulis, Wu, and Manifold worked on it 
for a small percentage of their time. 

0.5 Outline of the Technical Report 

This technical report is divided into two major parts: in Part I 

a numerical analysis of under integrated constraints is presented. Particular 



attention is focused on the Stokes problem with a constraint divergence 
u = 0 and on construction of an appropriate LBB condition for stability. 

Part II deals with underintegration and hourglass control. There 
projection methods, error estimates, and a large collection of numerical 


results are described. 



PART Is STABILITY OF SOME MIXED FINITE 


ELEMENT METHODS FOR STOKESIAN FLOWS 


1.1. Introduction 

In so-called primitive variable formulations of problems of flow of 
viscous, incompressible, Stokesian fluids, two fields appear as unknowns: 
the velocity field u and the pressure field p, the latter representing 
a Lagrange multiplier associated with the incompressibility constraint, 
div u = 0. Finite element methods based on such formulations were first 
introduced over a decade ago [42]. Since the mid-1970s, interest in these 
methods was rekindled by the appearance of several new techniques which 
provided for very efficient calculation of the element pressures. These 
included mixed methods which employ pressure approximations which are 
discontinuous at interelement boundaries as well as the closely related 
mixed-type methods which employ an exterior penalty approximation of the 
-incompressibility condition and reduced integration of the penalty terms. 

All of these methods have the attractive feature that the discontinuous 
element pressures can be eliminated element by element , reducing the problem 
to one only involving velocities. Upon determining velocities, element 
pressures can then be evaluated through a simple post-processing operation. 

Methods of this type were developed and discussed by several authors, 
and we mention in particular the works of Malkus [30,31], Hughes [23] » Malkus 
and Hughes j^3 b] , Reddy [jii} , Bercovier Q0» Engleman and Sani£l5]» and the 
references therein. In 1980, however, mathematical analyses indicated that 
some of the more popular discontinuous-pressure/mixed-methods might be numer- 
ically unstable (34 ,35-40, 4Q . It was discovered that while certain of these 
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methods perform well in problems with smooth solutions for which regular 
uniform meshes are employed, serious oscillations in the pressure approximation 
can occur when the data or the mesh pattern are mildly irregular, and these 
oscillations increase in amplitude as the mesh is refined. 

Oden, Kikuchi, and Song [41] attributed the deficiency of these unstable 

methods to their failure to satisfy a key stability criterion which they 

referred to as the "LBB-condition," making reference to the work of 

Ladyszhenskaya [29] on existence theorems of viscous flow problems and of 

Babuska [1 ] and Brezzi [8 ] on the approximation of elliptic problems with 

constraints. The discrete LBB-condition of Oden, Kikuchi, and Song is 

basically the requirement that the discrete approximation B* of the 

h 

transpose B* of the constraint operator B - div be bounded below as a 
linear operator mapping the space of approximate pressures onto the dual 
of the space of approximate velocities. For example, one form of this 
condition is that there exist an > 0 such that for all q^€ 


“hllOhM 


L (ft)/ker B* 
h 


sup 


( v dlv V 
IMi 


Related conditions for mixed finite elements were discussed by Fortin [L8] 
and Girault and Raviart [22 J. The possibility of unstable pressure approxi- 
mations is signalled by the existence of a parameter 

which depends upon the mesh size h. Indeed, the fact that a mesh-dependent 
“h corresponds to methods with "spurious pressure modes" is supported 
by the theoretical and numerical results of Oden et al [41 ] and by extensive 
numerical experiments of Malkus [32,]. Equally important, the behavior of 


* Definitions of terms displayed here are given in Section 1.3. 
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as a function of h governs the asymptotic rate of convergence of 
such mixed methods-. 

An important question that has arisen from these considerations is 
whether or not stable mixed methods exist which converge at optimal rates in 
the energy-and L -norms. The present paper is directed at resolving this 
question for a restricted class of problems by estimating the stability 
parameter in the corresponding discrete LBB-condition. 

The methods of proof of the LBB-condition basically fall into two 
categories depending on whether or not ker B* = ker B* or ker B*{£ ker B*, 

B* = — gradient + boundary conditions and B* is it finite element 

approximation. In the former case, a general method of proof can be constructed 

which is inspired by the work of Girault and Raviart [22], ^ndwhl.ch will be 

.discussed in the first section. We shall concentrate next on the latter 

case, and present another general constructive technique for estimating 

a, for uniform meshes which makes use of a discrete Poincare-type inequality, 
h 

These two methods of proof will be presented and used to establish the 

LBB-condition for two elements. In the first category of stable mixed 

* * 

method with discontinuous pressure for which ker B^ * ker B , we shall 
analyze in detail the Q^/P^-element (biquadratic velocity/ linear discontin- 
uous pressure) , and prove that it does in fact satisfy the LBB-condition 
with independent of h . 

As far as the second category of method (for which ker B^ (j£ker B^) 
is concerned, we shall prove that the I8/P^-element (eight node isopara- 
metric velocity/linear discontinuous pressure) satisfies the LBB-condition 

with ct of order h , and therefore appears to be unstable. However, 

’_-.h 

certain ways to stabilize this element are suggested and their implementa- 
tion in codes have lead to stable and accurate solutions. 



Also falling into this second category is the Q^/P^-element (bilinear 
velocity/piecewise constant pressure) • This element will be briefly discussed. 

Finally the various values of the LBB-constant and the rate of convergence 
expected from the most used rectangular elements and using discontinuous 
pressure will be summarized. 


1.2 Statement of the Problem 


Let ft denote an open bounded region of IR with boundary 3ft. 

We consider the two-dimensional Stokes problem an ft, which involves finding 
a velocity field u = (u^,u^) and a pressure field p such that 


\ 


-vAu + Vp = f in ft 
div u = 0 in ft 
u = 0 on 3ft 


( 2 . 1 ) 




where V is the viscoscity of the fluid, (v = const. > 0), and f is the 

2 

body force, assumed to be a prescribed vector field with components f^^ L (ft) 

We recast (2.1) in a weaker variational framework by introducing the 
spaces 

19 9 

( 2 . 2 ) 


v = (hJ(^)) 2 , Q = L 2 (ft) 


and the forms 


a:V x v -*■ ]R , f :V -*■ 3R 


> 


a(u,v) = v(u,v) , f(v) = \ (f.,v .) 

~ ~ x i=l 1 


(2.3) 


for all u, v€-V, where and (*»*) are inner products on V and Q, 


respectively, and are given by 
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\ 


(v,w) 


vwdx ; v,w€ 0 
ft 


2 3u . 3'v , 


= I ’ JT> ; u ’ v ^ V 

1 i,j=l 9x j 9x j 


(2.4) 


/ 


The partial derivatives in (2.4 ) 2 are interpreted in a distributional sense. 
We proceed by considering the problem of finding (u,p)€v x q such that 


a(u,v) - (p,div v) = f(v) Vv€v 


(q,div u) - 0 Vq € Q 


\ 


(2.5) 


It is easily verified that any solution of (2.1) satisfies (2.5); any 
solution of (2.5) satisfies equations of the form (2.1) in a distributional 
sense. Under the conditions stated, it is also known that (2.5) possesses 
a solution (u,p), with u uniquely determined by each choice of f and p 
unique up to an arbitrary constant. 

Problem (2.1) can also be interpreted as the characterization of a 
saddle point of the functional 


L:V x q -> 3R 

L(v,q) = |a(v,v) - f (v) - (q,div v) ( 2 . 6 ) 

with q clearly a Lagrange multiplier associated with the constraint, 
div v = 0 in Q. 

We remark the saddle point problem for the functional L(*,») of 
( 2 . 6), can be pre-conditioned by introducing the perturbed Lagrangian 

I £ : V x Q 4 i 

L £ (v,q) = L(v,q) -^|(q,q) 


(2-7) 
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for all q£Q, which represents a regularization of L(*,*) with respect 
to the multipliers q. For each e > 0, saddle points (u £ ,p £ ) of L(*,*) 
are characterized by 


a(u ,v) - (p^.div v) = f(v) Vv€V 


(eP £ + div u £ ,q) =0 Vq€ Q 

Upon solving the last equation in (2.8) for p £ , we obtain 

p = - — div u in Q 
e ~£ 


( 2 . 8 ) 


(2.9) 


The forms a(*,*) and f(*) are continuous and a(*,*) is V-elliptic. 

In addition, Ladyszhenskaya [29] has shown that a constant a > 0 exists 
such that 

cell q|| 2 < sup (q ^ lv 

L (ft) /3R v £ V II || 

M y If i 

Vq€ Q (2.10) 

where ||v|| = ~\/(v,v) 1 . Under these conditions, the sequence {(u ,p )} 

~ £ £ e > 0 

of solutions converge strongly in VxQ/H to the saddle point (u-,p) 
of the functional L(*,*) in (2.6). 

Finally, it is interesting to note that when (2.9) is introduced into 

the first equation in (2.8), one obtains 

a (u ,v) + — (div u .div v) = f (v) v £ V (2.11) 

which is equivalent to an exterior penalty formulation of the constraint, 


div u = 0. 



1.3 Finite Element Approximations 


We shall outline briefly features of certain finite element approxi- 
mations of (2.5) or (2.8). We confine our attention to cases in which ft is 
rectangular or is the union of rectangles and, for simplicity, to uniform 
meshes of rectangular elements of maximum length h. For a family of 
such meshes with E = E(h) elements, we introduce the discrete (finite- 
dimensional) spaces. 


vh - 


(v hl' v h2 ) 


v hi 6c (ft). 


v hi !f! « W : - 0 on 3!2 , 

e 

i < e < E, i = 1,2} (3.1) 

Q h - {q h €L 2 (® |, h | n 6P r (f! e )i 

e 

i < e<E, r>0) (3.2) 

Here Q (ft ) is the space of tensor products of complete polynomials in 
R G 

x, and x- of degree < k defined on finite element ft and P (ft ) 
is the space of complete polynomials of degree r defined on ft £ . The 
elements Q^/P^ and Q^/Pq clearly correspond to the values (2.1) and 
(1.0) of the parameters (k,r). 

In addition to the spaces V* 1 , we shall also consider cases in which 
V h is constructed using 18-elements: 

18 = eight -node isoparametric elements 
This element is also referred to as a serendipity element [46l . We also 
consider composite elements which employ both and 18— subelements. 



Clearly, for every h , 


V h C V and Q h C Q 


(3.3) 


The finite element approximation of the formulation (2.5) consists of seek- 
ing C V* 1 and p^ Q* 1 such that 


‘ <p h’ dlv V ‘ 

(q h ,div u^ = 0 


V v e V* 
~n 


¥<| h € q 


(3.4) 


while the approximation of (2.8) is of the form 


a(uf,v v ) - (p^,div v h ) = f(y h ) 


¥ V 6 v 1- 

-'ll 


rh’lh 7 vr h’ ~h 

(£p h + div Uh» q h* = 0 ^ q h € Q 


(3.5) 


The solvability of (3.4) depends upon a compatibility condition between 

the spaces V^ 1 and Q^ 1 which resembles (2.10) and which we record below. 

E £ 

Likewise, while (3.5) is uniquely solvable for u^ and for any £ > 0 

£ £ 

(under the stated conditions on a(*,*) ) , the behavior of and p^ 

as £ or h tend to zero also depends upon more delicate features of the 


app roxima t ion . 

* 

Let B, and B, denote the discrete operators, 
h h 

B h : vh "*■ Qh ; B h : qh ' - vh ' 


t q h’ B h^h ] = <V B h q h> E <V div ?h> 


Vq h € Q h , Vy h 6 V h 


(3.6) 


where [•»•] and , *^> denote duality pairings on Q' x Q (Q = Q' = L 2 (ft)) 

and V T x V respectively (i.e., B, and B* are the discrete approximations of 

n n 

div and -grad plus boundary conditions defined by (*,♦))• Then, the 
discrete LBB -condition for problems (3.4) and (3.5) is as follows: 

There exists a number a, > 0 such that 


q h^ L 2 (Q)/ker B* < s ^h 

h “ Zh €v 


( v div V 

Halli 


(3.7) 


for all q, 
n 



The behavior of as h tends to zero and the structure of ker B* 

governs the stability of these types of mixed methods. In particular, let 


E^(u,p) denote the distance function 


E^u.p) = inf || u - + inf || P - qjl 


f. v h 


q h € 


Q h 


L 2 (n) 


(3.8) 


defined on V x Q, § = (q € Q | J^qdx = 0}. Then one can show (see Oden 

and Kikuchi [40]) that if (u,p) is the solution of (2.5) and (u^,P^) 

h °b 

the solution of (3.5) in V x Q , 

II u - J^llj. £ C(1 + a" 1 ) (E^u.p) + e) 


II P - 




(3.9) 


where C is a generic constant independent of u, p, e, and h. 

The remainder of this part is devoted to the study of (3.7) and 
estimations of the stability parameter a, for different approximation spaces 
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V h and Q h ((3.1) , (3.2)). As noted in the Introduction, we will focus 
our study on the following approximations: 

1) Q^/P^ elements [biquadratic velocities, piecewise linear 

pressure^ 

18/P^ elements Ieight-node isoparametric elements for velocities, 
piecewise linear pressures] 

Composite elements [elements consisting of two or more of 
the abov^ 

q / Pq elements [bilinear velocities, piecewise constant pressures] 


2 ) 


3) 


4 ) 


Again, we note that in all of the cases we study, we shall assume that 


is a rectangle (or a union of 
rectangles) discretized by a uniform 
mesh of rectangular finite elements; 
V h CV, with V = (Hg(ft)) 2 and 
Q h d Q = L 2 (fl). 




/ 


( 3 . 10 ) 


The principal results concerning ker and the LBB constant oi^ are stated 
in the following theorems. 

Theorem J. Let conditions (3.10) hold and let the discrete spaces v h and 

& * 

Q h be constructed using -elements. Then ker B h = ker B and the 

stability parameter in the discrete LBB-condition (3.7) is a positive 

constant independent of h(a. • 0(1)) . 1 1 


Theorem II . 
defined by 
meter cl 


Let conditions (3.10) hold and suppose that V h and Q * 1 are 
18/P ^-elements. Then dim ker B^ * 3 end the stability para- 
in the discrete LBB-condition (3.7) depends linearly on h •* 


□ 


■ 0(h) 


are 


Theorem III. Let conditions (3.10) hold and suppose that V h and Q h 

2 

defined using composite 18/?^ - Q /P ^-elements of the type shown in Fig. 1 

£ 

Then dim ker = 1 and the stability parameter appearing in the dis 

crete LBB-condition is a positive constant independent of h ; (a = 0(1)). 

h 

Theorem IV. Under the assumptions of Theorem II 3 if V h and Q h are 

’k 

defined using Q^/P - elements > then dim ker « 2 and - 0(h) . 


1.4 The LBB-Condition For Q^/P^ Elements 

In this section, we describe a general method for establishing the 

k k 

LBB-condition when ker B^ and ker B coincide. This procedure will 

then be used to prove Theorem I for Q^/P^-elements, and can also be used 

for Theorem III. The method is embodied in the following four steps. 

I. Let q, be an arbitrary element in Q* 1 . Construct a vector 
h 

u^ £ v* 1 such that 


where || • |L - || • || 2 

L (SI) 


(q h ,div u h ) - II q h || 0 


llujli iCilqJlo 

and C is a constant. Then 


(4.1) 


sup 


<q h> divv h ) (l h .Jivu h ) 


» Fllijl 


V'i 


N ShNi 


h 0 


so that 


“h 


1/C . 


To construct such a u, , we continue as follows. 

~n 


* It suffices to define q^ only to within an arbitrary constant or to 
demand that all q^ be such that (l,q^) " 0 * 



II. For each q, C Q , q, f constant, it can be shown (Ladyszhens- 

h n 

kaya (_29J) that a £ V can be found such that 


div v q = q h in n and Hy^l^ q^ll Q 

Let w, denote the V-orthogonal projection of v, onto V : 
~n 

( r h ” Sh-^1 " 0 VY h € vh 


(A. 2) 


( 4 ^ 3 ) 


Then 


"hl'l i ll",Hl — C 1 II I' 0 


III. Set 


~ Yh ~h 


( 4 . 4 ) 


We attempt to construct a u^ with the desired property (4.1) by demanding 
that 

E f 

(4.5) 


JE. 

(q , div(e - e^)) = I q h div(e - e h )dx = 0 

e=l ■'fl 


where 


e, = u, - w 
~h ~h ~h 


(4.6) 


Then it is clear that 

( V div v h ) - II qj* - (q h ,dlv u h ) 

which is (4.1)^, and it remains only to verify that ( 4 . 1)2 ^olds. Assuming 
that this is possible, we see that the original problem reduces to one of 
constructing a u^ such that (4.5) holds. 
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IV. 


To satisfy (4.5), it is sufficient to require that 

h )dx 

+ I n*(e - e )dx 

-'an h ~ ~ ~ h 

e 

= 0 


'ft 


q, div(e - 


?h )dx 


-\ l 

J ft 


Vq, • ( e - e 
h ~ ~ 


for each finite element ft , n being a unit outward normal to 9ft . In 

e ~ e 

many finite element meshes, each ft^ is the image of a fixed Biaster element 
ft under an invertible affine map , 


F 

e 


ft ft^ , 


Fx = x = Tx + b 
e- - ~e 


(4.7) 


T being a 2 x 2 matrix and b a translation vector* Then it is suf- 
-e -e 

ficient to construct u, such that 
» 

q • (e - e, ) dx - o A qn • (e - e, )dfi =0 (4.8) 

J ft ~ ~h Jan ~ ~ ~ h 

where q = q,oF \ e = eoF \ etc. 
h e ~ ~ e 

Remark: This procedure is next used for the Q^/P^-element . But in 

the case where ft is partitioned into I8/P^ elements, except one Q^/P^ 
element, we can show that 


^ & 
ker B, = ker B 
n 

For this mesh, the construction II, III, IV can theoretically be made, 
but the essential estimate of (4.1) cannot be obtained. This remark suggests 
the introduction of the composite element described in Theorem III. 

For the the two discrete finite-dimensional spaces V* 1 

and Q* 1 are defined as 
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yh ‘ { ih ' <v hl’ v h2 ) l v hi € c °® ; 


v hJn e «2 <n . ) > v hi ' 0 on 


an, I < e < E, i = 1,2} 


Q h - \ € L 2 («) | . qJ Q £ P 1 (n e ) I < e C E} 

e 


and then using the definition 

* r h f 

ker B = {q, C_ Q such that q, div v, dx = 0 
h h . h ~h ~ 

for all v, €. V h } (A. 9) 

n 

a simple calculation reveals that 

ker B* = ker B* = 3R (4.10) 

n 

Then it suffices to construct a u, such that (4.8) holds for the 

~h 

/\ 

master element ft shown in Fig .2 and to then show that u^ satisfies 
(4.1) £ • We use the notation indicated in the figure; the integral appearing 

~ /V 

on the left side of (4,8 ) is denoted I , and we seek u with u^ £ Q^Cft). 
Observe that the shape functions associated with the indicated nodes are 
of the form 

1|> 5 = (1 - x 2 )(l - y 2 ) , = %(x 2 - l)y(l - y) 

^ 23 = %x(l + x) (1 - y 2 ) , $ 34 = h(l - x 2 )y(l + y) 

= %x(x - 1) (1 - y 2 ), etc. 
q € P^(^) is of the form 


and that each 
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q = q 0 + q i 5 + q 2 y w±th ^ q = (q-^^) 

where q a » ot = 0,1,2 are real numbers. A simple calculation reveals that 


I = q 0 [ 

(a - g) 

/V . A 

•n ds 



°-'9(2 





+ 

/A 

-v ^ e hl 
(2 nx 

e^) dxdy + 

4 

A/ A, /V \ A jA-1 

x(e, - e) *n dsj 
9Q ~ h ~ ~ 

(4.11) 

+ q 2 C- 

4 

- (e. 7 - 

(2 

e^) dxdy + 

y(e, - e) *n dsj 
9(2 ~ h ~ ~ 



It is clear that we can make 1=0 by choosing e^ (equivalently, 
choosing a u^ ) such that the following five conditions hold: 

(i) e^a 1 ) = 0 r 1 <i<4 

(ii) e^Ca 1 ^) = 0, 1 <_ i j <_ 4 

A /v 

where T is the unit vector tangent to 9(2 
e^*n ds = | e*n ds, 1 < i £ j <4 

e^dxdy + | xe^*n ds = -| e^dxdy + o xe»n ds 

9(2 (2 9(2 

e, „dxdy + I ye -n ds = -f e dxdy + o A ye-n ds 
nz J a ~h ~ 1 Jao — 


(Hi) 

(iv) 

(v) 


(2 


(2 9(2 (2 

This set of conditions must determine the 18 independent components of 


% «h± £ V a)) 


Conditions (i) and (ii) make 12 of the 18 degrees of freedom of e, 
zero. We are left with six coefficients: 


\l * *hl ( * 5) ®5 + *hl C ‘ 23) ®23 + \l (al4> *14 

hi ' 4 h2 (a5) < , 5 + *h2 (al2) 'i'l2 + *h2 (a34) 'i i 34 
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of these coefxicients 3X6 iiuiQ6di3t6ly determined from (iii) by 
a direct integration: 


- / 12, 

a h2 (a > 


a 


A » A 

e^ds ; 


a , 23, 

e hi (a > 


' ?} A dS 


~ , 34, 

e h2 3 > 


3 

4 


A m A 

e 2 d S ; 


A , 14, 
e hl (a > 


a 

■ - t{ A" 


Thus, it remains only to determine e h± (a 5 ), i = 1,2, using the last 
two conditions, (iv) and (v) * But a direct calculation leads to the 
pair of equalities 


a a 

i|e h i(a 5 ) - J^dSdy + J ^dS 


S ,*V* - l\ 


a 

2 


16a , 5, 

-S e h2 (a > 


Ae dxdy - if e n ds - f yg ds 

fi - 1 1 ^ J 2 1 


4 


A • A 

e 2 ds 


d 

AA . A 

ye, ds 

J 9 x 


Kence, conditions (i)-(v) determine a vector g for which 1=0 

-h * 

required. We easily verify that e €v h ,e, I = e. oF . 

-h -h ~h e 

e 


as 
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It remains to be proven that the vector u, = e, + w, satisfies 
* ~n ~h -n 

(A.l)^* note t * lat ^ i s sufficient to prove that 

II ~tJI 1 — C i II ~ II 1 (4 - 12) 

because 

fl-hiii- H 5 .H 1 - |,u h - " iA — °xH z, - aA 

ivAA + AA> 

so, since || <_ || v 1^ 

na.Hi « (l + jyii^i < cm , h n 0 

To establish (4.12)* we note that for the master element** 

§h = iwA » II Shl'i ” - c * II I ! h ( ^i ) II l 2 (4 - 13) 

where {b^} are t ^ ie 9 nodes of the element and || | • j| | denotes the 

2 

euclidean norm in 1R . Using the fact that 

ij^e-n dsi < cnaio.as < ccn«n5, fi + ml#' 


lj a i*e-4 ds | < C|| S|l 0 


an 


, etc. 


and the previously computed nodal values of e^ obtained via steps (i)- 
(v) above, we can verify that |] | ± ) U | 2 ± c|| e|| 2 jj + II ®ll^ q- Thus, 


a constant C exists such that 


A *hHi, nic{ listing HI »llj f a> v 


(4.14) 


* Here and elsewhere in this paper, C denotes a generic constant independent 
of h and does not necessarily have the same value throughout. 



We next transform this result so that it applies to a typical element 

of the mesh and sum over all elements to obtain 
e 

IISh' l i,8- cth ’ 2|l ? ll o + (4 ' 15) 

In this last calculation, we used the affine map of (4.7), the fact 

that 11 | T e II l C^li, 11 1 T e ~ 1 11 I c 2 h_1 > and standard relations between 
|| ||| 1 _g and || ;|| 1>a . 

We shall next verify that 


II e|l Q <_ Ch |1 e\\ 1 


(4.16) 


We will then arrive at (4.12) via (4.15) and thereby complete the proof of 
the theorem. 

To prove (4.16), we employ a duality argument of Girault and Raviart 
[223. Note that 


(e ,v) 

II e i II o a = s ^ p 2 — ITV * 1 " x » : 

* V £L 2 (fi) Ilv|| 0#n 


(4.17) 


Let g be in L (ft) and <J> be the solution to the Dirichlet problem 


• s 


<t> E I an ' 0 


Then 


(4.18) 


♦ Cn 2 (ft) PiiiJ(a) and II *gll 2 , n i c ll silo, n (4,19) 


The variational formulation for the problem (4.18) is: 

( V v, i,a ’ <8,v) o,n 

It is permissible to take v = e^^ so that (<fr .e^)^ ^ = (g, e^Q 
But e i is orthogonal to V h ; hence, (v^.e^ ft “ 0 » V v h 


Q 


0 



It follows that <e i>g ) 0>!! - (e ± ,* g - V in V v h «■ V a " d 

l( ®**i>o,n l - 11 'A, n 11 * 8 ' v J l i,a V- h C'' h • 

Choosing v, = 4*^ to be the interpolant in V* 1 of d> , we have 
n g g 

ll»8-*8 ll l.^ C hH*ll2, n i Cb » SM 0 , n 


Hence, 


„ „ 

— — £ Ch!|e || and by (5.7) , 

8 ' 


0, (2 


n-iio.oi-«*n -n x , a 

This completes the proof of the theorem. 


1.5 The LBB Condition For I8/P^ Elements 

This section is devoted to the proof of Theorem II and in particular 

* 

to obtaining the kernel of the operator B and an estimate of the LBB 

h 

constant a 
h 

vh ■ { ’h ■ (v hl> v h2 )|v hi£ c ° (fl) ! 

v hi l 0 £ q 2 (! V' v hi ■ 0 on 

e 

3ft, 1 £ e £ E, i = 1,2} 

= Q 2 C^ e ) " {x 2 y 2 ,X €. 1R, (x,y) £. ft g } 

Q h - tq h Q|q h l n £ 1 £ e £ E} 


where is the subspace of used to define the serendipity element. 

/s. 

We shall work with a master element ft . Each element q £ Q has 
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three degrees of freedom , a = 0,1,2, and each element and can be 
chosen such that 

q[ e = q Q + q x x + q 2 y ; Vq = (q^ q 2 ) 

A 

When referred to the element , these degrees of freedom will be noted 

by q a , a = 0,1,2 : 

$lfi " ^0 + V* + V ; 3 = ( V V 

Then 

q 0 = q 0 ’ q l = hq l * 3nd q 2 = hq 2 * 

* 

Lemma 5.2. Under the conditions of Theorem II, dim ker = 3 .. 

Moreover 3 ker B, = span {l, x, » X,J » where X-, » X„ aTe discontinuous 
n 1 2 12 

functions of the type shown in Fig. 3. ... 

Proof . It suffices to confine our attention to the collection of four 

A 

reference elements i = 1,2, 3, 4, shown in Pig* 4a. 

We wish to characterize all q^ £. ker B^ ; i.e., (q^> div v^) ■ 0 

\/v h V* 1 . We begin by choosing in V* 1 such that =0 at all 

14 14 

nodes except a where v^(a ) = 1 • Then, 

in n i : v hl = -4x(l + x) (1 - y) , q h = qj + q*(x + j) + q^Cy - j) 

in : v h4 = -4x(l + x) (1 + y) , q h = q£ + q*(x + j) + q^(y + j) 

0 

where q , a = 0,1,2, are the degrees of freedom (coefficients) of q, 
u h 

A 

in subelement , e = 1,2, 3, 4. Then we find that 



q, div v, dx 
n au 


0 
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and this implies that 


q l = " q l 


14 


Similarly, choosing v h2 (a ) = 1 with = 0 at other nodes, leads to 
the conclusion that q* - q^. We continue this procedure at nodes a 12 . 


a 23 , and a° H to eliminate 7 of the 12 degrees-of-freedom of q h - Collecting 
these results, we are left with the 5 coefficients indicated in Fig. 4b. 

We next choose - 0 at all nodes except that Vj^O) = 1* We 
find that 


34 


’ 

q.div v, dxdy = 0 => q 0 = q~ 

/s a /» h ~h l l 

n-U n 2 U fi 3 U n 4 

A similar calculation with = ^ y^ e ^ s *lx - ^x* 

Collecting all of the results, we are left with three independent 

coefficients, <! 0 ><Jx’ q 2 and these define the q h -pattern indicated m 
Fig. 1. Reciprocally, a linear combination q^ of l»Xj and X 2 satisfies 
(q h ,div v h ) = 0 for all y h in V h . A typical member x> of ker B £ is 
indicated in Fig. 3 ; X 2 is obtained b y rotating the x-, y-axes 90-degrees. o 


Proof of Theorem III . We now return to the completion of proof of 
Theorem II. On each element , we evaluate the product q h div y h using 

16 degrees of freedom of y h (eight for each component v ± and v 2 > and 
the 3 degrees of freedom of q^ (the coefficients q 0 ’ q l’ and q 2 de ^ ined 


earlier) . We get 



q h div v„ d* = J 

a 

e 

+ \ v^a 3 ) [+q o +2q 1+ q 2 ] + \ v^a 3 ) l+VV 2 ^ 1 
+ | v^a 4 ) [-q o + 2 q r q 2 ] + f v^a 4 ) [+^-^+ 2 ^] 

+ j [~v 1 (a 12 )^ - v 2 (a 12 )q Q ] 

. 2 r , / 23, * f 23. ✓. n 

+ 3 [ +v 1 ( a )Q 0 - v 2 < a )q 2 ^ 

2 r , 34. /s , , 34. ^ , 

+ j L-v 1 (a )q 1 + v 2 (a )q Q ] 

2 r , 14. /v ( 14. /s 

+ 3 t_V l (a )q Q -v 2 (a )q 2 ]J 

E r 

(q h » div Vj,) “ I j q h div v h dx 
S=1 ft 


{ ? V* 1 ) + ? V al) m 0 +«i+j« 2 ] 

+ ? - 1 (a 2 >[+q 0 +2$ -g 2 ] + i v 2 (a 2 )[-8 -q +2$ ] 


This summation over the elements can be replaced by a summation over 
the nodes as shown in Fig. 5. Three types of nodes can be distinguished as 
shown in the Fig. 6 . The indices for the pressure with respect to each node 
is also shown. Then, using the numbering scheme shown, we have 


2 h _ 1 (q h , div v^) 


1 r / / V r /si - /si ,/sl , ^2 , /s2 a2 , a3 , «a3 , a 3 /1^4-9n^ n^l 

= 6 i ( v i( e i ) I- £ lo +2q l +q 2 +q o +q r q 2 +q o +2q l +q 2" q o +2q l q 2 ] 


+ v 2 (d x ) 

+ I £ ( v l (e II )[ -V^ 1 + v 2 (e XI ,[ -^ + ’o 1 ) 

6 II 

+ 3 ^ ( v l (e III )I -^o 1 + v 2 (e III >t 'V^2 1 ) 


'III 



If we choose: 


•• / ^ \ /T / ^ ^ I A ^ I j ^ I i A 2 * aa2 A 2 | A 3 | a A 3 | A 3 A^ • a A ^ A ^ \ 

W - 6(-q o +2q l +q 2 +q o +2q 1 -q 2 +q o +2q l +q 2 -q o +2q 1 -q 2 ) 
v 2 ( Sl ) = 6(-q^q x +2q 2 -qJ-^+2q 2 +q\q\2q^-qJ+2q^) 

v 2< e ii> ' I HfiJ) 

v 2 (e H ) ' I 

v l (e UI ) ■ 2 ^VV 

v 2 ^IIl' 3 7 ^2 _q 2^ 

and u = 0 at the nodes on the boundary , 


(5.1; 


then 


2h 1 (q h » div y h ) 


N 


L( v i (n)2 + v 2 (n)2 ) - c 1 1 ~h 1 1 1 


n=l 


(5.2) 


where the summation is over all N nodes. 


Now it will be shown that the choice (5.1) implies 


HaJll - c llq hll 0 /k e r B* 

n 


(5.3) 


Then (5.2) and (5.3) complete the proof. D 


2 

With the expression chosen in (5.1) we can reorder ll v hlli an d> using 


the numbering scheme of Fig. 5, obtain 
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1 1 Zh 1 1 1 1 c I ( (-%+2qi+q2^^q^?^q^q^q^^ 2 qJ-^ )2 

+ (-qJ + ^i + 2 t 2 "^o" 2 ^i + 2 ^ 2 + % + ^i + 25 2 + ^o“^i + 2 ^ 2 )2 

I ( r^“ 4~n^ '\^ m 4- / , /v 3 \ 2 » /^l,^2v 2 ✓ /-v 3 , /v4 N 2 

+ vq 1 + q 1 J + (q i +q i ) + (q 2 +q 2 ) + ^ q 2 +q 2 ) 

4- n ^ ^ ^ i /a2 a1 \ 2 r/v3 a4\2 a.^4 a1 \ 2 *\ 

+ (q o -q o ) + (q o -q o ) + (q^) + (q^) ) 


(5.4) 


Here we find a quadratic form, whose kernel is precisely the kernel 
of B* defined in the Lemma 5.1. 

Then, it can be written: * 


II || 2 s _ r f tf3- ~2,2 , ( f2 .3,2 ..3 .4,2 ,.4 .1,2 

" 5 i"l 1 C I { ( V q o } + (q o“ q o ) + (q o -*c? + ‘VV 

i ^ 

X 1 A A A n A - 1 A < 


,.l .2,2 .2.3,2 ,.3 .4,2 ,.4 .1,2 

+ ^-q^ + (q^) + (qj^) + (q^i) 


(5.5) 


■ y . 3 , . 2 , 2 . ...2 . 3,2 / a 3 , . 4 \ 2 , . 4 A 

+ (q 2 + q 2 ) + (q 2 -q 2 ) + (q 2 +q 2 ) + (q 2 - q 


i> 2 } 


The passage from ( 5.4) to ( 5.6) comes from the fact that both 
quadratic forms in bracket in these expressions provide the same kernel 
and therefore define two equivalent semi-norms on 

V = (q* ; i = 1, 4 ; j = 0, 2) 

C J 

Now if we pay our attention to the quadratic form on the first line 

2 

of ( 5.5), it can be interpreted as the L -norm of the gradiant of a 
piecewise bilinear function <J>q , defined by q* at the corresponding 
centroid of each element. This function <{)q belongs to and, as 

proved in TEMAM [45], there exists a constant C such that 
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ll 7 » 0 ll 2 i c||* 0 || 2 

L /IR 


(5.6) 


This procedure can also be applied in the construction of functions 

<j>^ and 4>2 from the ' s and ' s. Finally we obtain the inequality 

(5.3) from (5.5 ) and (5.6 ) noticing that the summation of the three 
2 

squared L / IR norms of <{>_. (j ■ 0, 1, 2) precisely corresponds to 



2 


0/ker B 


* 

h 


1.6 The LBB Conditions For Q^/P Elements 

The general proof of Theorem III can also be used in the analysis of 

the Q 1 /P q element (bilinear velocities, constant pressures). For this element 

& 

we maintain that the kernel of B h consists of checkerboard nodes which are 
characterized by alternating values a and b in each neighbor element. In 
this case 

* 

dim ker B =2 
h 

Using the same notation as in the previous section, we can define for 
each element and each node their integer component J and K ; the element 
e = 1 in the corner (resp. e = 2) corresponds to J = K » 1 (resp. J = 2, 

K * 1) • (See Fig. 7.) Using the elements e satisfying J + K e even and 
constructing a piecewise bilinear continuous function defined by q at 
the centroid of these elements, we can apply the inequality (5.6) and obtain 


N 

I 

i=l 

J+K even 


,°3 


*3 2 

+ 


i djE 

i=l 1 
J+K odd 


*4.2 

q i ) 


Ch 


E 

l 

e=l 

J+K even 


° 2 


We obtain a similar inequality considering the element e such that 


J+K is odd, and by addition, we arrive at 
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N 

l 

i=l 


,»1 3,2 . ,o2 »4,2 

^ q i “ q i^ + ^ q i " q i^ 


Ch 


f ° 2 

\ q e 
e=l 


( 6 . 1 ) 


2 4 " 

The improvement h instead of h in the estimate of Oden, Kikuchi, and 
Song [41/] allows us to obtain an LBB constant = 0(h) for this element. 

1.7 Summary of Some Stability Results 

A mathematical analysis of the discrete Babuska-Brezzi condition (3.7) 
has been made by Oden and Kikuchi Bs] , Oden, Kikuchi and Song Bq] , and Oden 
and Jacquotte ^ or sever al finite elements for a model two-dimensional 

Stokes' problem on a uniform mesh. We shall summarize these results here 
which pertain to the behavior of the "LBB-constant” and the stability of 

the pressure calculations. We use the notations 

= space of complete piecewise polynomials of degree k 
over an element 

Q k = space of tensor products of complete polynomials of 
degree k 

18 = the eight-node isoparametric element 

Results are summarized in Table 1. In this table, examples 1, 2, and 7 
"lock" at small values of the penalty parameter e . This means that for a 
given mesh size h , e cannot be taken arbitrarily small, as noted earlier. 
Of course, for an acceptable e for reasonable mesh sizes, e is so large 
that the constraint of incompressibility is not adequately satisfied. Hence 
these elements should generally be avoided. Elements 2, 4, 5, 8, 11, and 
14 are unstable since = 0(h) . Remarkably, these instabilities frequently 
are not observed on uniform meshes when the solution is very smooth. Mild 
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irregularities in the solution or small perturbations in the mesh may, how- 
ever, produce violent oscillations in computed pressures the magnitudes of 
which increase without bound as h tends to zero. In many cases, however, 
these oscillations disappear upon "filtering” the pressure solutions (i.e. 
upon averaging the pressures over one or more elements) ; In the case of 
elements 2 and 14 it has been proven mathematically 20 that certain filter- 
ing schemes will produce a stable and convergent method. However, it is not 
known if filtering can be used to stabilize and salvage the remaining unstable 
elements. 

Elements 6 and 10 lead to stable and convergent schemes and are quite 
robust in the sense that they are insensitive to singularities in the solution. 
However, they are not too accurate and converge at a suboptimal rate. 

Element 9 is clearly the superior of any listed: it is unconditionally 

stable, it provides both velocity and pressure approximations which converge 
at the optimal rate, and 

ker V - ker V 
n 

Element 13 is somewhat of a novelty. While element 5 yields unstable 
pressure approximations, Oden and Jacquotte Qf] have shown that a composite 
of three e ^ ements ( no * 9) and one I8/P^ element (no. 5) is stable. 

The behavior of elements 11 and 12, marked with an asterisk, is only 
conjectured here and has not been rigorously proven. 

Extensions of these results to three-dimensional elements are straight- 
forward. 

1.8 Numerical Examples 

The results of several numerical experiments are described which are 
designed to verify the theoretical results with regard to the Qj/P^, 
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and the composite elements described earlier. We also investigate numeri- 
cally the effects of a pressure filtering operation. 

As a first example, we consider an L-shaped domain ft partitioned into 
64 square subdomains, as shown in Fig. 8. The fluid is subjected to a 
constant body force f = (0,-100). We take V » 333, and the penalty 
parameter e = 10 We will be interested in the computed hydrostatic 
pressure across the section AA f defined by: y * 0.80 . Each subdomain 

corresponds to a finite element; the velocity on each element is interpolated 
at 8 or 9 nodes and the pressure by its value at 3 points. Thus, various 
choices of how to handle the ninth node lead to meshes with 18/P^, Q 2 /P^, 
or Composite/P^elements. We will be interested in three cases involving 
these elements: 

Mesh 1: All the elements are C^/P^ elements 

Mesh 2: All the elements are I8/P^ elements 

Mesh 3: Adding 16 centroid nodes, we obtain 16 composite 

elements as shown in Fig. 9. 

The results reported here were obtained using the FIDAP code for problems 
of incompressible viscous flow 0-4] „ 

Figures 10 and 11 show the comparison between the results obtained with 
the Q 2 /P 1 element (Mesh 1) and those obtained with meshes 2 and 3. Fig. 10 
illustrates the major difference between the C^/P^ the I8/P^ element; 

the former involves a pressure which seems to be smoothly distributed along 
the section AA 1 while the latter yields a pressure with severe oscillations. 
We note, however, that the values of the pressure obtained at the centroid of 
each element are close to the values obtained with the Q 2 /P^ e ^ ement » which 
suggest that this unstable solution can be stabilized by a filtering opera- 
tion which effectively uses these averaged values of pressure. 
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It is also remarked that the oscillations seem to come from the 

spurious modes in ker B*. The smoothing device may be equivalent to an 

a posteriori elimination of these spurious modes. However, it turns out 

that these spurious modes do not solely come from ker B* : Figure 11 

contains results obtained by adding one node in the elements in the corner 

(point in Fig. 8). For this mesh, ker B* = 1R but the results still 

exhibit pressure oscillations. However, for this mesh, the solution seems 

to be much smoother than in the 18-case. 

Finally, the composite elements lead to a quite smooth solution as 

indicated in Fig. 12, which is close to the solution obtained with 9-node 

2 

elements, except that for this element h = 0.25, while for the (^/P^ 

2 

element h was equal to 0.0625. 

We also note that when the body force f derives from a potential: 
f = _Vv , then the unique solution for the Stokes Problem is 

( u = 0 
P ■ -v 

In this example f = (0,-100) and v = lOOy = -p. The numerical results 
obtained by these different methods are summarized in Table II. We 
observe that the error in the filtered pressures is around 33 percent 
greater than that of the Q 2 /P elements for this particular example. 

As a second example, we consider a Dirichlet Stokes Problem, which is 
designed for numerical verification of the convergence theory for three 
schemes considered: 

a)_ A uniform mesh of e ^- ements * 



b) 


A uniform mesh of 


I8/P^ elements. 

c) A uniform mesh of 18/P^^ elements with pressure averaging. 

We consider the unit square domain partitioned into- square subdomain, and 
the following body forces f = (f^fj) are applied: 

f = -4y + 12x 2 + 24xy + 12y 2 - 24x 3 - 48x 2 y - 72xy 2 

-8y 3 + 12x 4 + 48x 3 y + 72x 2 y 2 + 48xy 3 - 24x 4 y - 48x 2 y 3 
-2 (x - x Q ) + a(x) 

f - 4x - 12x 2 - 24xy - 12y 2 + 8x 3 + 72x 2 y + 48xy 2 + 24y 3 
-12y 4 - 48xy 3 - 72x 2 y 2 - 48x 3 y + 24xy 4 + 48x 3 y 2 

where a(x) = -1 if 0 £ x £ x^, a(x) = 1 if (Xq < x £ 1. Then 
(u,p) is defined by 

/ 

= X 2 (l - x) 2 (2y - 6y 2 + 4y 3 ) 

< 

u 2 = (~2x + 6x 2 - 4x 3 ) y 2 (1 - y) 2 

2 

p = x Q - X - (jc - x Q ) if 0 £ x £ x Q 

< 

2 

P - x - x Q - (x - x Q ) if x Q < x £ 1 


u | =0 

~ r 

s div u = 0 in Q 


u = (u^.Uj) ; 


and 


(u,p) satisfies: 


-Au + p = f in Q 
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As before, we construct a plot of the computed pressures across a 

section of the domain. Figure 13 shows the results obtained by partitioning 

the domain in 64 square subdomains. For this mesh, h is equal to 

1/8. The computations are made with Q2 -on 18-elements. Whereas the C^- 

solution seems to be stable, clearly the 18-solution shows oscillations 

around the exact solution. However, it is noted that both solutions 

coincide at the centroid of the elements and this again suggests that 

the "smoothed 18-solution," obtained using only the pressure at the centroid, 

2 

is stable, and may converge at a rate of 0(h ). 

Finally, Fig. 1!4 confirms this suspicion showing the computed rate of 

2 

convergence is precisely 0(h ) for the pressure for the (^-element, and 
for the smoothed 18-element. However, it is also observed once ae;ain that the 
Q^/P^-pressures are considerably more accurate than the filtered I8/P^- 
pressures for all mesh sizes considered. 

With the results from these examples we can conclude that 

2 

• The C^/P^ element is stable and the optimal L -rate of convergence 

2 

of the pressures of 0(h) is attained. 

• The I8/P^ element yields unstable pressure approximations, but 
these can apparently be stabilized considering only the values at 
the centroids. 

• Spurious oscillations (checkerboarding) can also appear when 

ker B* = 1R 
n 

• Filtering the pressures in the 18/P^-element by using only the 

centroidal value leads to a pressure approximation which may 
2 2 

converge in L at a rate of 0(h); however, the accuracy of the 
filtered scheme is quite inferior to that of the (^/P^-elements. 



These computed results underline once again the critical role played 
by the LBB-condition in studying the stability of finite element schemes 
by reduced integration. These and other results we have computed also 
indicate that the estimates obtained in Section 4 for the discrete LBB- 
constant are sharp. Indeed, the theoretical result that the use of a. 

composite element of the type employed here leads to a stable pressure 
field, while not of great practical value, is fully confirmed by the numer- 
ical results. This suggests again that these calculated estimates of ct 

h 

are a good indication of the actual numerical performances of these methods. 



PART II: ANALYSIS OF INSTABILITIES IN 


UNDERINTEGRATED FINITE ELEMENT METHODS 
2.1 Introduction 

For many years, a special type of numerical instability has been 
observed in finite difference approximations of flow fields, which has 
been referred to as "hourg las sing", "keystoning", or "chickenwiring". 
These graphic terms refer to geometrical patterns which appear in com- 
puted flow fields (e.g. velocities) and which emerge as spurious oscil- 
lations superimposed on an otherwise smooth field, the spurious oscilla 
tions often taking a zig-zag form which resembles an hourglass or a 
chickenwire mesh. These spurious modes can be amplified upon refining 
the mesh, and to control such numerical instabilities, various schemes 
for incorporating "hourglass viscosity" or "hourglass damping" have 
been proposed by some authors. 

It is now known that hourglass modes can arise from an incomplete 
(or poor) approximation of the kernel of the operators in the momentum 
equations in flow or solid mechanics problems (or, more generally, of 
the principal part of the operator in the governing differential equa- 
tion of a given boundary-value problem) . For example, in addition to 
the rigid body motions residing in the kernel of the standard operators 
appearing in the equilibrium (momentum) equations of solid and fluid 
mechanics, one finds hourglass modes in various crude discrete models 
of these operators. 

In recent years, the occurrence of hourglass instabilities in 
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under integrated finite element approximations has been observed . In 
the implementation of most finite element methods, integrals defining 
stiffness matrices are evaluated using numerical quadrature schemes. 

To improve computational efficiency, the practice of undBT^tnt&g^cztion 
is often employed, by which is meant the use of a quadrature rule of 
an order lower than that required to integrate polynomial integrands 
exactly. This can produce rank-deficient stiffness matrices or, equi- 
valently, an expanded kernel of the equilibrium operation which contains 
spurious hourglass modes, and the result is again a numerically un- 
stable scheme. 

In order to overcome this difficulty, artificial stiffness or vis- 
cosity methods, or other stabilization methods have been proposed by 
several authors (e.g. [2—5, 17, 28]). These methods involve computing 
an under integrated matrix, and then adding a stabilization matrix which 
effectively eliminates the hourglass modes. They turn out to be fairly 
general and have been used for a long time in numerous codes. Whereas 
all of these methods based on intuitive feeling give good numerical 
results, their mathematical study remains often non-existent. 

The most interesting challenge is to solve the problem using only 
the crude rank-deficient underintegrated stiffness matrix, the solution 
is obtained up to within an arbitrary spurious mode, and then to elimi- 
nate these modes from the solution in a post-processing operation. 

Unfortunately, even when the stiffness matrix is rank-sufficient, 
similar oscillations are observed when under integration is used. In 
that case, the process of the excitations of modes similar to the hour- 
glass modes is not completely understood and these modes have never been 



mathematically studied. 

In this report, we would like to give precise mathematical justi- 
fications and answers to the questions previously mentioned. The 
next section (Section 2.2) is devoted to the proof that the Stabiliza- 
tion Method is mathematically justified. Then, in Section 2.3 we pre- 
sent a method which involves solving an underintegrated and not well- 
posed problem, then in a-posteriori eliminating the unknown degree of 
freedom. The proof of the accuracy of the method is given in Section 
2.4, and its numerical aspects and results are described in Section 2.5. 
In Section 2.6, we examine the case where the spurious oscillations 
cannot be predicted from the rank-deficiency of the stiffness matrix and 
we analyze why these modes may be excited. Finally, we apply the pre- 
vious considerations to the Linear Elasticity Problem. 

It should be noted that the method and its results cannot be 
embedded In a classical elliptic theory: Strang f s ellipticity condition 

[44] is here violated and this non-elliptic method cannot be studied by 
the classical theory of finite element methods and numerical integration 
[10, 11]. We also refer to Girault [20, 21] for his approach to the 
same kind of problem, non-elliptic because of the use of partially 
under integrated stiffness matrix, but where hourglass modes did not 
appear. 

2.2 A-Priori Hourglass Control 

2.2.1 Intr oduct ion . This section is devoted to giving a mathema- 
tical support to several methods consisting of adding a stabilization 
matrix to the underintegrated matrix. For clarity, we shall confine our 
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attention to a simple model problem. Let ft be a regular domain in 
]R with boundary 3ft and consider the model Neumann problem. 


<v 

Find u = u(x,y) 

such that 




-Au = 

f 

> 

on ft 

(2.1) 


3u m 

3n 

0 

on 3ft J 

i 


2 

where f is an L (ft) -function satisfying 

[ f dxdy =0 (2.2) 

Jft 

The questions of the existence and uniqueness of solutions to (2.1) 
(which are well-known) are taken up in Part II. 

We shall first consider a finite element approximation of (2.1) 
constructed using Q^— elements, i.e., four— node quadrilateral elements 
over which bilinear shape functions are used. Most of our notations 
and results are reproductions of those of Flanagan [17] and Belytschko 
[2, 4, 5]-- Then we will attempt to extend our results to the Q_- 
elements (nine-node, biquadratic elements) and will indicate in which 
ways they differ from Belytschko's [3]. 

The construction of finite element approximations of (2.1) involves 
the calculation of the stiffness matrix for a typical finite ele- 

ment ft , which is given by the formula, 
e 

K = VN C * Vn dxdy (2.3) 

~ e Jft ~ 

where N is a vector representing the bilinear or biquadratic shape 

functions in each element ft , 1 < e < E . 

e — — 



When (^-(respectively Q^-) elements are used to discretize the 
domain £2 , is a 4x4 matrix (resp. 9x9) and the N f s contain four 

bilinear (resp. nine biquadratic) shape functions. We will distinguish 
exact-, full-, and under- integrations . The full integration is obtained 
using the number of Gauss integration points necessary to obtain the 
exact integration on regular square elements: 4 (resp. 9) points in our 

study. The underintegration will involve the Gauss rule of lower order: 
1 (resp. 4) points. The stiffness matrix associated with a rule involv- 
ing k points will be denoted , k = 1,4,9 . 

Several authors [2, 4, 5, 17, 28] proposed to add to the underinte- 
grated stiffness matrix a stabilization matrix which exhibits several 
special properties. In this section, we will prove that these proper- 
ties are indeed satisfied and that the exact stiffness matrix K can 

e 

be computed by this method. This will be accomplished by first carry- 
ing out the integration (2.3) exactly. 

We first introduce some notations. Suppose that element £2^ is 
defined by the coordinates of its nodes (x*, y 1 ) , 1 < I < p, p = 4 or 
9 . We introduce the isoparametric mapping from a master element 



to £2 such that 
e 


x = 


l X 1 N T (S,n) 

1=1 

> 

p 1 
l 7 N 

1=1 J 


(2.4) 


y 
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where Nj , 1 I p , are the shape functions for the quadrilateral 
element on the master element. The node numbering convention is shown 
in Figure 15. 

The stiffness matrix K is evaluated in (2.3) using the mapping 

~e 

/v 

(2.4) from £2 to £2 : 

e 


K 

~e 


J £2 



VN dxdy 


J £2 


A J 

J| VN 


"d£ 

H 

r&i 

dx 


a. 

X 

l 


VN d£dn 




> 


(2.5) 


is the Jacobian matrix of the mapping from to £2^ , 

J is the inverse of its determinant, and where the gradients of the 
shape functions are derived with respect to the master element coordi- 
nates (£,n) • These matrices can be computed using (2.4) : 


where 


d£ 

dx 



dx _ 

_ dC J 




dx 

dx 


d£ 


T 

dN 

T 

v 

dN 




d? 

T 

dN 

T 

dN 


dr) 

y 

dn 

T 

dN 

T 

d!J 


dn 

-y 

it 

T 

du 

T 

Y 

dN 


dn 

A 

d ? 


( 2 . 6 ) 


(2.7) 


( 2 . 8 ) 
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VN 




(2.9) 


where J is the Jacobian of the mapping 


J 


det 


dx 

dZ 


( 2 . 10 ) 


Finally, we obtain the expression. 


K 

~e 


. T T . T. T' 

Axx A Ayy A 

+ 


'ft 


y T Ax 


T. 

y Ax 


|d?dn 


( 2 . 11 ) 


where A is the antisymmetric matrix 

A _ d|J dN T dN dlj T 

~ = dn d? ‘ d? dn 


( 2 . 12 ) 


the Jacobian J can be expressed as y Ax . 

A study of expressed as in (2.11) and the properties, of the 

matrix A will then enable us to study the effect of the underintegra- 
tion of the stiffness matrix. We will first concentrate on the 4 node- 
element and derive the exact expression of the stabilization matrix. 
Then we will discuss what form this matrix may take for the 9-node ele- 
ment. 


2.2.2. The stabilization matrix for the bilinear element . For the 
bilinear element, the shape function vector can be written as 


N 




+ Sn h 


(2.13) 
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where 


s 


8 - 
t T = 

t T = 
h T = 


a. -i, -i. i) 
a, i, -i, -i) 

(i, i, i, i) 

(i, -i, i, -i) 


(2.14) 


then the explicit form of the (4x4) matrix A is 


= i(s's T - ss ' T ) +|(sh T - hs T ) + £(hs ,T - s'h T ) 


A » (w w 
4 - - 


4 




(2.15) 


for (£,n) ■ (0, 0) , we obtain which satisfies 


?V ■ r!| r •„ 5 ' ^ y 24 x i3 + 'nV <2 - 16) 

• £=n=o 

which is merely the area of the element , noted |ft e | . 

At this point, we can precisely see the matrix resulting from a 
1-point rule; this underintegrated matrix denoted by is given by 


. T a T . T T 

„(D . ^oH^o 

K e ' “1^1“ + “l^T 


(2.17) 


Also, if we note B = (b^, b^) the discrete approximations of 
the gradient VN evaluated at the integration point we can remark that 


Si ■ -ttol 


S 2 ■ So: 


(2.18) 


and therefore (2.17) takes the usual form 


5e 1>= ( SlS? + S 2 S 2 > 


(2.19) 
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the rank deficiency of can now be verified. 

Indeed, from (2.14) and (2.15) we can simply note that 


Jo?" ° 


“ 0 


( 2 . 20 ) 


and then 


K h = 0 


K t - 0 


( 2 . 21 ) 


Therefore, if we consider H and T the global hourglass and transla- 
tion, and the assembled under integrated stiffness matrix, we have 


•H = Z *H 

e 


also K 


Z K (1) »h 
~e <» 
e 


( 2 . 22 ) 


and this proves the rank deficiency of K v . . 

Note that this ”±1" pattern is independent of the regularity of the 
mesh and that H will take alternating values +1 and -1 at neighbor 
nodes as shown in Fig. 16. 

Our goal will now be to calculate a matrix K Sta ^ such that, if 

~e ’ 

added to , we obtain the exact stiffness matrix given by 

(2.11). This expression does not seem easy to integrate, but the image 
of certain vectors mapped by this matrix can be easily computed using 
orthogonality relations previously obtained (2.20) and the fact that 


= b^y = |nj 
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T T 


(2.23) 


we can obtain: 


K t - 0 


K x = b 
~e~ -1 


(2.24) 


5e? ' *2 


Equation (2*24) is not sufficient to compute , because it gives only 

9 out of the 10 coefficients of K (4 x 4, symmetric) * It is enough 

~e 

to know X t K X , where t, x, y, and X form a set of independent 
~ ~e ~ ~ ~ - 

vectors* That is the case for X * h because 


det(x,y,t,h) = 4A ^ 0 


provided the element is not singular* Then the knowledge of h K^h 
and the relations define uniquely . If we set 


h K h = 16 e 
- 


(2.25) 


then is given by 


(1) — T 

= K v; + e YY 


K = K N 
~e ~e 


(2.26) 


whereas, again, given by (2.25), e is a scalar, and 


T t 

Y = h - y b 

X 5 n Ei nT 2j 

A 1 


(2.27) 


While it is difficult to express £ nicely in function of x and 


y , its exact value can be written 
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- 1 
4 


[s(s T x) - n(g ,T x)] 2 -t- [c(s T y) - n(s’ T y)] 2 

lO.I + ^y 4 3 X 12 + y 12 X 34 } + n(y 32 X 14 +y l4 X 23 ) 


d^dn 

(2*28) 


We observe that for parallelogram elements, the denominator is constant 
and its value is the area of the domain . In this case 


24 K 


, 2 . 2 . 2 . 2 . 
x 13 + x 24 + y 13 + y 24 ) 


(2.29) 


or for rectangular elements 

l 2 + l 2 

F- — * 

12 £ l 
x y 


(2.30) 


where £ and £ are the lengths of the sides of the rectangular 
x y 

element. Also note that for such parallelogram elements y reduces to 
h . 


The expression (2.26) is often used to a-priori eliminate spurious 
modes for the kernel of K but the determination of £ remains a prob- 
lem. The choice e * 0 leads to the underintegrated matrix and to the 
method to be studied in the next section. On the other hand, a cheaper 
way than the full integration of the whole matrix would be to fully 
integrate £ given by (2.28). This method would lead again to the 
full integration and is cheaper because it needs only one 4x4 integra- 
tion by element instead of .10. A more common use is to take for £ a 
simple value independent of the geometry of the element, which is often 
the value obtained for a square 1/6 or sometimes any arbitrary constant, 
as used in [4, 5]. 



2.2.3. The stabilization matrix for the biquadratic element. In 


this section, we will study the effect of the underintegration of the 
9x9 stiffness matrix obtained in (2.11) with nine-node elements when 
a 4- Gauss integration point rule is used. Whereas Belytschko and al. 
[3] intuitively obtain another "y.y 1 " stabilization, we prove that this 
decomposition is not even valid for regular meshes. We then propose a 
decomposition derived on regular mesh. 

(4) 

But first we exhibit the spurious modes out of . For the bi- 

quadratic. element, the shape function vector can be written as 

N = S | (2.31) 


where S 


is the 9x9 matrix 


S 



-2 -2 4 

2-2 4 

2 2 4 

-2 2 4 

0 4 -8 

—4 0 —8 

0 -4 -8 

4 0-8 

0 0 16. 


and 


f = [l, S, n, Cn, S 2 , n 2 , 5n 2 , 



(2.32) 


(2.33) 


The integration rule we are interested in involves four integration 

points (£ a , n a ) , a = 1, 4 • Associated to each of them, we note 

A and J the corresponding matrix A and Jacobian, The underinte- 

~a a r 

grated matrix is then 
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„< 4 ) i + &a\ 

K — L 

~e - J 

a=l a 


(2.34) 


and can also be written 

4 


r ( 4 >- I i (b?bf + b?bf) 

- e 0=1 J o - 1 - 1 - 2 - 2 


(2.35) 


where 


B a = 


“1 

, (XT 


, v T 

b -l 


(-V) 

. OtT 


T 

(A a x) 

-2 




(2.36) 


generalizes (2.18) to a 4- point rule. 

(4) 

The rank deficiency of K g can now be verified. Indeed if we 
call t and h the vectors defined by 


T 

t = 
h T = 


, 1 , 1 , 1 , 1 , 1 , 1 , 1 , l] 

[l, l, l, i, -l, -l, -l, -l, o] 


(2.37) 


we easily obtain: 


and 


T T 

t • N = t • S • £ 


h T • N = h T • S • £ = -4(£ 2 + n 2 - 12 C 2 ti 2 ) 


and then differentiating these expressions and using (2.12) we get: 


A ■ t = 0 


A-h - -8 £i( 1-12£ 2 ) || + ?(l-12n 2 ) ||] 


the second expression vanishes when the point (£,n) is one of the four 

A 

Gauss integration points of ft 
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( 5 °, ti®) = (± -7=-, ± -i— ) ; o - 1 , 4 
2V3 2V3 


therefore we have 


(2.38) 


A 

-a 


*0 ' 5 - 0 


a = 1, 4 


and consequently 

K (4) .t-K (4> -h-0 


~e 


(2.39) 


(2.40) 


which proves the rank deficiency of 


( 4 ) 


Once again we note that 
the pattern of h defined in (2.37) is independent of the geometry 
of the element and is therefore valid for a rectangular element as well 
as for an irregular element. 

The search for decomposition for this element cannot be completed 
in a manner as complete as it has been for the 4-node element. However, 
Belytschko and co-workers [3] have intuitively come up with a decompo- 
sition similar to (2.26) where y and e are 


1 h T 

4 * 


b a 


4 b a 


~1 

J 

1 . T 

- 4 * * 1 

112 

e , f 

(2.41) 

a 


a=l a 



y s\ 

h * *2j 


(2.42) 


This decomposition does in fact satisfy several properties also satis- 
fied by the exact matrix, as 


K . t = 0 
e 


K . x = 7 b“ 

~ e - £1 - 1 
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K 

~e 


4 

y - z 

a=l 



but for a simple square element , K g and its decomposition (2.26) does 
not coincide. Indeed, for this simple geometry, the calculation of 
(2.11) can be carried out explicitly and the polynomial in (£,n) ob- 
tained can be split into one part exactly integrated with 4 Gauss points, 
and another part of higher order that requires 9 , points. This calcula- 
tion leads to the decomposition: 


r (9) = k (4) 

XX XX 

r (9) = k ( 4) 
yy xx 


+ n 

e 

+ n 

e 


,1 T . 4 T. 

( 45 £7* ?7 + 135 £9 * & 

, 1 T 4 T. 

( 45 !8* & + I35 f 9 * 




> 


J 


<2.43) 


where 


K 

xx 




A 


T 

• y • y 


7 Ax 



d£dn 





> 


y 


(2.44) 


and , Sg , are the 7 t * 1 , 8^ and 9*"^ column vector of S (2.32). 

These vectors correspond to the higher order of £ (2.33) that cannot 

be exactly integrated by a 4-point rule. The form taken by the stabil- 

T 

ization matrix involves now three matrices (s^* , i = 7,8,9), is 

exact for a square element and cannot coincide with the decomposition 
found in [3]. Finally, we note that both decompositions were used in 
our a-posteriori control described in Section 2.7 on a regular mesh. 


* or also for a geometry for which the Jacobian is constant. 
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and optimal rates of convergence were only obtained with the decomposi- 
tion (2.43). 

2.2.4. The stabilization matrix for a general heat transfer equa- 
tion. In this paragraph we would like to give the stabilization matrix 
for a slightly more complicated operator. The case of the linear elas- 
ticity operator will be discussed later. 

Let us consider the case where the operator is defined by 

(2.45) 



A = 

T 

B C 

B 

where 

B - 


_9\T 

3y ; 

and 

C - 

f c n 




\ C 12 

C 22 1 

then the 

stiffness 

matrix associated 


* - [ 

T 

VN ■ 

1 C • VN dxdy 


(2.46) 


n 


The generalization of the stabilization decomposition when elements 
are used can then be written 


K_ = + e y • y T 




(2.47) 


where 


K U) 

~e 


a 


T 

B C B 


(2.48) 


e ■ C,, e + (c, _ + C-.)e + c _ 0 e 

11 xx 12 21 xy 22 yy 


(2.49) 
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and 


- 1 

e = ~r 
xx 4 


n 


IT T 2 

•j(£s.y - ns’* y) dEtfn 


yy 


xy 


1 

4 

!_ 

4 


~ ^r(?s T * x - r|S ,T * y) 2 d^dri 

n J ~ ~ 


> (2.50) 


/v - 7 (£s T • x - ns T * x)(£ s?y - tis ,T • y) d£dn J 

Q J 


The quantities Y, B and J are the ones previously defined. Expres- 
sions similar to those given in (2.29) and (2.30) can be used to 
simplify £ . 

For a regular geometry, and corresponding to (2.29) and (2.30) we 

have 


Sex " 24 (ft ) (y 13 + y 24 ) 

e 


1 , 2 , 2 x 
Vy“ 24(jy (X 13 + X 24 ) 


Sty " 24(ft g ) (x l3 y 13 + X 24 y 24 ) 


(2.51) 


As far as the 9-node element is concerned, the decomposition can 
be obtained only for regular elements* First we note that 


k < 9) . j-C4) 

- xy - xy 


(2.52) 


where the notations are similar to (2.43) and (2.44). Therefore, the 
decomposition can be written: 

K< 9 ) - v(4) . A r , l T 4 T. 

-e ? e Q e C ll ( 45 ?7 ?7 + 135 ?9 ?9* 

+ S e C 22 ( 45 ?8 ?8 + 135 ?9 ~9 ^ 


(2.53) 



2 >3 A-Posteriori Hourglass Control 


2.3.1. Introduction and preliminaries. The basic ideas are more 
easily understood when demonstrated for the same simple model problem. 
We still focus on the model Neumann problem Pq or its variational 
equivalent P . 

2 

Let ft be a regular (e.g. Lipschitz) domain in H with boundary 

2 

3ft and let f be a given L -function. Problem P^ is then, 

(Pq) Find u such that 


-Au = f in ft 


= 0 on an 
3n 


(3.1) 


where the data f satisfies the compatibility condition , 

f fdx = 0 (3.2) 

J n 

Later we shall put further restrictions on fi and on f (e.g. we 
will need f 6 H(ft)). The kernel of the governing operator A = (-A , 

g 

•s— ) in (3.1) is, of course, the space of constants. Thus, whenever 
on 

(3.2) holds, there exists a solution to (3.1) which is unique up to an 
arbitrary constant. 

To formulate a variational statement of problem P^ , we introduce 

* 

the spaces and inner products , 


2 

* The elements of V (and L (ft)/H) are cosets C v 3 such that u G [yl 
implies that u, v G H^(ft) (or L^(ft)) and v - u G H . Throughout this 
paper we frequently refer to functions v in V , meaning, of course, 
that v is a representative function in the coset . 
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V = H 1 ^)/ TL 


(u,v ) 1 

C f »g) Q 


an inner product on V 
Vu*Vvdx ;u,vGV 

a 

2 

an inner product on L (Q)/ M 




fgdx - 


meas 


J 0 £to i 


fdx I gdx 

a 


(3.3) 


Three remarks are in order: 

i) The norm || • ||q associated with the inner product 

2 

the canonical norm on the quotient space L (fi)/ H , 

||f || - inf || f + X|| 

x e ® l (fi) t 


ii) According to Temam [45 ], there exists a constant C Q , depend- 
ing only on ft , such that 

IMll l c 0 lMlo Vv 6 V ' (3.5) 

iii) For all f satisfying the compatibility condition (3.2) and 
any v G V , we have 


(f.v), 




fvdx <_ II f || Q ||v|| 0 

it l|fII 0 II- 111 


(3.6) 


With these relations now established, we consider the variational 


statement of Pq as problem P : 



57 


(P) Find u 6 V such that 

(u,v) 1 = (f,v) 0 VvGV (3.7) 

We can easily verify that any solution of Pq is a solution of P 
and, conversely, the solution of P satisfies the condition of P^ in 
at least, a distributional sense. Moreover, since the bilinear form 
(• > *)^ is continuous and coercive on V and since the linear form 

is continuous on V if (and only if) f satisfies (3.2), the 
following result is an immediate consequence of the Lax-Milgram Theorem: 
THEOREM V. Let f satisfy (3.2). 

Then there exists one and only one solution u 6 V to problem P 
and this solution depends continuously on the data f . □ 

We now consider a finite element approximation of the problem P. 
Let us now construct a finite element approximation of problem P. We 
begin by introducing a partition Q of f2 into E finite elements 
so that 

E 

n = U 

e=l 

We shall assume that ft is such that it can be partitioned in this 
fashion into four-node quadrilateral elements over which bilinear shape 
functions are defined. Thus, if 

Q (ft ) = space of bilinear functions defined on ft e 

we can introduce the finite-dimensional space 

V h = {v^ 1 6 C^(^) such that 6 ) , 1 fi. e E / II C. V 

e 6 (3.8) 



wherein, as usual, the label h is the mesh parameter (e.g. , 

h = max dia(Q )). The functions in V* 1 are continuous and are still 
l<e<E e 

defined up to an arbitrary constant. 

Our finite-element approximation of problem P is embodied in the 
discrete problem, 

(P, ) Find u* 1 G V* 1 such that 

" (3.9) 

(u h , v h > 1 - (f, v h > 0 Vv h e v h 

where, again, f satisfies condition (3.2). 

In analogy with Theorem V, we have: 

THEOREM VI. Let f satisfy (3.2). Then there is one and only 
one solution u* 1 to 'problem P^ in V* 1 and this solution depends 
continuously on the data f . Q 

In examining the convergence of such finite element approximations, 
we shall confine our attention throughout this study to regular mesh 
refinements. In such cases, we have the a priori asymptotic error 
estimates, 

|[ u-u 11 II x -0(h)-, II u-u h || 0 = 0(h 2 ) (3.10) 

2.3.2 The underintegrated problem . We now focus our attention 
on finite element approximations of problem P in which incomplete 
quadratures are used to evaluate the bilinear form (*,*)^ . To simpli- 
fy this study, we shall now introduce some additional assumptions: 
i) Q is the unit square. 
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n = (0,1) x (o.i) 

ii) The finite elements are the squares, 

= (iii — ) x ( -i) 

“ij ^ N * N ; K N ’ N ; 


1 < i , j < N 


l<i,j<N 


iii) The data f is L -integrable; e.g. 


f « L (Si) 


(3.11) 


In this case, we take 


h - i , dim V h = (N+l) 2 - 1 = 0(h' 2 ) 
N 


(3.12) 


^ 2 

In P, we can replace f by f , its L -projection on V is 
n 


defined by 


(f\ v h ) 0 = (f, v h ) 0 


Vv h 6 v h 


(3.13) 


For further use, we note that the projection satisfies 


II £ h ll 0 < l|t|l 0 


(3.14) 


and can be chosen such that 


f h dx = 0 


(3.15) 


Now we turn to the issue of numerical integration of the stiff- 
nesses. Let !(%•) denote a discrete inner product on C^(fi) defined 


by a numerical quadrature rule as follows: 
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E 

I 6 (f.g) = 1 I®(f,g) 

e_1 (3.16) 

I®(f.g) - 2 W®f(^)g(?, e ) 

e j-1 J ~ J 

0 g 

Here W_j are the quadrature weights and £ are the quadrature points 
for element e and G is the number of quadrature points used. 

Assuming that Gaussian quadrature is used, the choice G=4 (2x2 - Gauss 
rule) leads to an exact integration of the stiffnesses for each element: 


, h h. , h h. T 

(u , V ^ = I 4 (u , v ) = u K v 


(3.17) 


for any u* 1 , v* 1 6 V* 1 . Here K is the fully-integrated stiffness 
matrix and u and v are vectors of nodal degrees of freedom of u* 1 
and v* 1 , respectively. 

Instead of the correct bilinear form in (3.18), we wish to consider 

an under integrated approximation to (•,•) in which only one integra- 

★ 

tion point per element is used: 


, h h. _ x / h K T._(l) 

( u » v h X 1^ U » v ) = u K v 

\J h h _ „h 

V u , v 6 V 


(3.18) 


.( 1 ) 


Here K is the underintegrated stiffness matrix. The difference be- 
tween ( # , - )^ and (*,*)^ ^ (on V* 1 ) is denoted a'(*,») and the 


corresponding stiffness matrix is K' 


stab ** 


* Recall Section (2.2). 

Recall that K Sta ^= £' 
Y is given by~(2.27).‘ 


** Recall that K = £YY where £ =1/6 for a rectangular mesh and 
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, , h h. , h tu , h h. 
a ’ (u , v ) = (u , v ) 1 - (u , v ) ljh 


T Tr stab 
u K v 


w h h _ „h 
V u , v G V 


(3.19) 


The "underintegrated problem", 


(P. ) Find u G V such that 
n 


<“ h ■ A,h - < fh ■ v ’ h 6 vh 


(3.20) 


is, in general, meaningless. This problem, in general, has no solution 
except for the special case in which f* 1 is orthogonal to the one- 
dimensional space of hourglass modes. 


H = { H G V h |(H, v h ) L h = 0 V v h G V h ) 


(3.21) 


A way to overcome this difficulty is to note that the under integration 
of the righthand side also leads to a rank— deficient linear form 
(*>*), 


0,h 


(f h , H) = .0 , V f h g V h 

V H 9 h 


Note that if i satisfies (3,15) we also have 


< fh • l >0,h • 0 


(3.23) 


Therefore we now consider the underintegrated problem P^ : 


(P ) Find ^ 6 V* such that 
h 

, — h h N 

(u • V >l.h (f ’ V } 0,h 


V v h G V* 1 


(3.24) 



where 


v * 1 : V h /H 


We can now state and prove 

THEOREM VII : There exists one and only one solution u* to T . ' 

h 

Proof : This is an immediate consequence of the Lax-Milgram theo- 

rem. Since 

(u h , uh ) 1>h = 0 < ^ > u h = y 1 H + y 2 

we can consider (*,*)- . as a norm on V , It is therefore coercive 

l,h 

and continuous on V . As far as the continuity of the righthand side 

h h 

is concerned, a simple calculation shows that for any v in V we 


have 


« h . v h , 0 >h l< ||f"|l 0 ||v“|| 0 


Also for any constants y^ and y^ 


r h h 


(£", v" + Y 1 + y 2 B) 0>h | - |<f\ v h ) 0>h | 


therefore 


II v" + Y x + y 2 h || 0 W y 2 
il|f h H 0 l|v h + Y 2 H|l 1 Vy 2 

l^lU h H 0 Hv h |l 1>h 


Here we successively used (2.23), (2.22), (3,6) and the equivalence 


— h 

between the canonical norm of V and the norm || *|| 


l,h 


□ 


We have obtained a solution to the under integrated problem P 



This solution is unique in V , from a computational point of view it 
is defined up to within an arbitrary hourglass mode. We now need a 
projection to obtain a reasonable solution from any representative u* 1 
chosen. 


2.3.3. Projection of the underintegrated solution . In order to 

construct this projection, we remark that, since u is a solution of 

P, and since H 6 V* 1 , u* 1 satisfies 

h 


(u h , H) x = (f h , H) 0 


(3.26) 


We wish to extend u* 1 to all of V* 1 so that a new function u h 
G V* 1 is obtained which contains an hourglass mode and which also satis- 
fies (A. 8). Thus, if tt is an operator from V* 1 into V* 1 , we define 


u^ 1 = nil* 1 * 5* + XqH, X q G ]R 


(3.27) 


({i h ,H) 1 = (f\ H) q 


This latter requirement determines X Q uniquely as 


X. = 


° 


[(f h , H) q - (^ , H)J 


so that u* 1 is uniquely determined as the function 
-h -h , ( fh - H >0 n „ 

Ml II H I!? 


(3.29) 


It is instructive to consider a geometrical interpretation of our 
projection defined in (4.9). Note that the "component" of the fully 
integrated solution u^ orthogonal (in V ) to H is (u ,H)^ = (f »H)q, 
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as indicated in Fig. 17. The solutions u, of P. constitute 

n h 

the vectors generating a line "parallel to" the space H in the figure. 

projection u is then the vector defined by the orthogonal projec- 
tion of u onto this line. Indeed, by construction. 


(u h - U h , H) 1 - 0 

At this point, we have established the following procedure for 
processing an under integrated finite element approximation of problem P. 

i) Compute the underintegrated bilinear and linear forms (• •) 

h 1,h 

and (f . •) 

v ^ 0,h 

ii) Solve problem P^ for u 
Hi) Compute (u* 1 , H> 1 

iv) Construct the enhanced solution u* 1 using (3.30). 

Thus, this procedure involves the computation of an underintegrated 
solution u^ to a reduced problem P^ and its enrichment via a post- 

i. 

processing operation to obtain a new approximation u . We shall now 
show that these post— processed solutions u^ 1 converge to the exact 
solution u of problem P as the mesh is refined, and, remarkably, 
these approximations converge at precisely the same rate as the fully- 
integrated solution! 

Indeed we have: 

THEOREM VIII : Let u , u* 1 and u* 1 be the solutions of P, 

P, and P , let f be in L 2 (fi) and satisfy (3.2). Let u* 1 be ob- 



defined in (3.30). Then we have the 


tained by the projection of u 
following error estimates for s = 0 and 1 

||u h -S h || e <C x h 2 * * - s ||£|| 0 (3.31) 

and 

ll u - ah |l s h 2 " S Ilf II 0 (3.32) 

The next section will prove this theorem, 

2,3, Convergence of the A-Posteriori Control 

2.4,1 Introduction , This section is devoted to the proof of 
Theorem VIII. The method of proof relies on the tensor properties of 
the bilinear element and of the Gauss integration rules. The problems 
and will be explicitly solved using an orthonormal basis of 

eigenvectors of (*,•)-,» (•,•'), . and (•,•)«. . Then we note that 

i- -L ^ ri (J y n 

2 2 

for a regular domain and mesh, f 6 L (ft) implies u € H (ft") and that 
II u - » h ll 1 < Chllf Il 0 (4.1) 

Likewise, the Aubin-Nitsche method provides also 

II u - u h || 0 < c 'h 2 ||f !| 0 (4.2) 

By the triangle inequality, 

II u - 5 h H 1 < Ch ||f|| 0 + ||u h -u h || 1 (4.3) 

with a similar estimate in the ]| ‘HQ-norm . 

Thus, it suffices to estimate the relative error 
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h 

e 


~h 

u 


h 

- u 


(4.4) 


2 1 

The L - and H -‘norms of this relative error will be explicitly calcu- 
lated and estimated. 


2.4.2. Some one-dimensional results . For reasons to be made 
clear in the next subsection, it is convenient to review briefly some 
results on one-dimensional piecewise-linear approximations on a uniform 
mesh for SI * (0,1) . Our aim here is to establish concrete relation- 
ships between various bilinear forms (*,*) n u > (•,•)«, , and (*,*) 1 

u,n x x, n 

on spaces of piecewise-linear functions. 

Let D(k, a) and I denote the N+l-order matrices 



"k 

a 

• • 

• 

o 

o 

1 


~ 1 

0 

• 

0 

0 ~ 


a 

2k 

» « 

o 

o 

• 


0 

2 

* 

0 

0 

p(k,a>= 


• 

• • 

• • 

il f - 


• • • 

• 

• 



0 

0 

• • 

• 2k a 


0 

0 

• 

2 

0 


0 

0 

• • 

• a k 


0 

0 

• 

0 

1 





" 






(4. 

(i.e. 

I' 

= (D(l, 0) . 

Then, for 

o^O 

, one can show 

that 





det D(k, a) 

=* (-a) N+1 det D(- 

k 

a » 

-1) 








- (-a) N+1 det D(- 

h 

a 




(4. 

where 



def 










D(k) 

= D(k, 

-1) 






(4. 


The values of k for which det D(k) vanishes are 


k - cos , 0 <_ i £ N 


(4.8) 
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and the corresponding vectors (D(k i )v i = 0) are 


v, = (cos -5-} , 0 < j < N 
~i N — — 


(4.9) 


The significance of the above matrices is that in one dimension, 

12 2 
the discrete H (0,1)- , L (0,1)- and underintegrated L (0,1) -norms, 

on the space of piecewise linear C^-functions on a uniform mesh 

of N elements on (0,1) , 


« {v* 1 6 C°(0,1) j v* 1 is linear on [eh, (e+l)h], e=0,...,N-l} 


are associated with the matrices 


(4.10) 


A Q = | D(l,%) , A x = ^ D(l, -1) and k Q ^ = | D(l,l) (4.11) 


respectively. In other words. 


hii 2 

v = v A v 
11 s ~ ~s~ 


s = 0, 1 , (0,h) 


(4.12) 


where v is the vector of nodal values of v 

By using (4.6) through (4.8), one can verify that the numbers a. 
and $ i which render A Q - c^Aq and ^ - 6^ singular are 


3(1 + cos -^) 
1 2(2 + cos ^) 


. iTt 

6 1 - cos T 

1 h 2 2 + cos 

N 


In particular, let (j>* = <{>*(x) , x G [o,l] denote the piecewise 
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linear functions associated with the vectors : 


<t> 1 ( jh) = cos , 0 £ i, j < N 


span 


(4.15) 


Then, 


t h / 11 i i. 

(v > * Vh “l (v ■ ♦ >0 

(v h , ^ - 8 ± Cv h , 4> 1 ) 0 


h _ „h 
v G V 


(4.16) 


(4.17) 


Notice that the base functions <}> are orthogonal for each of the 
scalar products under consideration. 

The following remarks are in order: 

i) The denominators in (4.13) and (4.14) are non-zero, 

ii) For i = N, = 0 and the corresponding eigenfunction is the 

one-dimensional hourglass mode: 

( 1 , - 1 , 1 , - 1 , ...) 

iii) For i = 0, g ■. 0 and the corresponding eigenfunction is 

constant. Then we have the condition (v^,l)^ * 0 as expected. 

2.4.3. Discrete norms for two dimensional meshes . The extension 
of the above results to two-dimensional rectangular meshes is straight- 
forward. Since the bilinear basis functions for V* 1 are tensor pro- 
ducts of piecewise linear functions of one variable, we can define 

^(x.y) = 4> ± (x) (y) 


0 < i,j < N 


(4.18) 
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Further, let us normalize these basis functions so that 


IU 1J II 0 - 1 

We can then establish the following: 

Lemma 4.1. For v h G V h , we have 

, h A i -U , h « ij * 

<v • ♦ >0,h ■ “i “j (v > * >0 

<v\ ♦ 1J ) 1 - (B ± + Bj) (v\ * iJ ) 0 

* 13) i,h ‘ <a j 6 i + “iV * 1J >o 

Moreover 3 if arbitrary v* 1 G V* 1 is expressed in the form 3 


h v , ij ^ 

v = Z v 4> 

0<i,j<N 1J • 


, h i ij v 
v ij * < v > * >0 


> 


then 


H E 

0<i, j<N J 


||v h ||J = Z (B + 6 )v^ 
0<i, j <N J J 


Proof: First note that 


, h ij. 

(v » <f> ^ 0,h 


T , h i j . 

I(V <J> (x) <t> (y)) 

a . a . [ ■ v h ^(x^y) dxdy 

1 


- Cj (V h . *«) 0 


(A. 19) 

(4.20) 

(4.21) 


(4.22) 


(4.23) 

(4.24) 



70 


We also have 




9yk ,i ' . . , j . . 

<)> (x) <r(y) 

+ (y) <t )i (x) dxdy 


= v^^CxJ^Cy) dxdy 

+ & I v h( )> i (x)<^(y) dxdy 
^ ^ fi 

= (& ± + Sj) ( v h , 4> ij ) o 

Finally 


(vh * ♦ 1 J> i,h = vfr * 1,<|,j + It ♦ i * 3 ' ) 


= S i a j (v h ,<J ) :i:3 ) 0 + B j a i (v h t <{) i:i ) 0 


The norms (4.23) and (4.24) are then directly obtained ^ 

In analogy with our remarks on the one-dimensional case, we observe 
that for i = j = N, = H , the two dimensional hourglass mode. Then 

, h . 24 , h 

(v , H) 1 = — (v , H) q (4.25) 

h 

and 

<A H) o,h = 0 <4-26) 

Also, for i = j = 0 , <j>^ = 1 and the equilibrium condition (3.2) can 
be written 



(4.27) 
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with 


f ij - « h - * iJ ’o 


(A. 28) 


2.4.4. Explicit resolution of and (P^ + tt) . With the above 
results in hand, let us now return to the fully-integrated finite- 
element approximate problem given in (3.9). The solution u* 1 to 

that problem can be written 


u 


0<i,j<N 


"ij* 


ij 


u ij = (uh * ^>0 


and since for [jv* 1 = in (3.9)], 


(4.29) 


(u h , <{> ij ) 1 - (B ± + Bj) (u h , <j) ij ) 0 
“ (fh - * 13 >o • £ li 


we have 


u 


ij B ± + Bj ij 


f.. ; (i.j) f (o,o) 


(4.30) 


Using constructions similar to those in (4.29) for the fully- 

integrated problem, we easily verify that the solution u* 1 to the 

underintegrated problem P, is representable in the form, 

h 


-h 
u = 


E u..* ij 

(i,j)*(N,N) 1J 
(i,j)^(0,0) 


(4.31) 


u . . = 


a jJi 


1 J a i Bj + ij 


£,. ; (i.j) ¥ (0,0) and (N,N) (4.32) 


with 



The cases (i»j) — (N,N) and (i,j) — (0,0) correspond to the arbi- 
trary hourglass mode and arbitrary constant, respectively. 

The projected approximation u h defined by u* 1 = ttu* 1 is con- 

h h 

structed so that projections of u and u coincide; i.e. 


"ij = U ij * <°>°> and (N » N) 


U N,N " U N,N 


(A. 33) 


2.4.5. Proof of theorem VIII. Since the error function e h = u h - 


u h is in V h , we use (4.29) and (4.31) to obtain 


l e . .<p 
(i,j)?f(N,N) 1J 
(i,j)? £ (0,0) 


(4.34) 


where 


e ij - (e • ♦ >0 


= <u\ - (u\ <|> ij ) 0 


U lj " U ij 


Thus, from (4.30) and (4.33), 


a i a j 


ij 


a. 6. + a.B. 

i J J i 


B, + B, 1 f ij 


(4.35) 


Then, using (4.13) and (4.14), e^ can be written as 


e , . = h k f 
ij *^ij *ij 


(4.36) 


where 



(4.37) 


K ij “ K < cos Y * cos ) 


and 


K(x,y) 


1 (1+x) (1+y) 

4 (1-bc) (l-y)+(l+y) (1-x) 


1 (2-hQ(2+y) 

6 (2+x) (l-y)+(2+y) (1-x) 


(4.38) 


On the square S = ^-1,+lJ x C~l> + l] /{(-1,-1) , (1,1)} , K(*,*) is 
bounded and there exists a positive constant R such that 


|K(x,y) | < K V(x,y) 6 S 


(4.39) 


Therefore we have 


l K ijl i K V(x,j) ^ (0,0) and (N,N) 


(4.40) 


and we can obtain using (3.13) and (4.23) 


h„2 


h E 

(i»j)^(0,0) 

(i,j)^(N,N) 


2 2 
K ij f ij 


Ih 4 K 2 !| f h || 2 < h 4 K 2 1| f || 2 


Also, after calculation and use of (4.24), (4.14) and (3.14), we have 


h II 2 

e "l 


h E 

(i,j)^(0,0) 
( i, j)^(N,N) 


+ V K u 


£12 h 2 K 2 ||f 11$ 


2 . 5 . Implementation and Numerical Results of the A-Posteriori Control 
For the Laplace Equation , 


In this section we first would like to indicate how the a-posteriori 



control method is implemented, and how its time efficiency compares to 
the a-priori method. Then several numerical results will be given, 
illustrating the accuracy of the method and confirming the results 
obtained in the previous sections. 

2.5.1. Implementation of the a-posteriori method . First let us 
indicate that from a mathematical point of view the problem P* 1 is 
well-posed but computationally, the matrix obtained from this formula- 
tion is singular and the dimension of its kernel is 2. Consequently, we 
must pick two nodes, fix them a value, and solve. The first value fixes 
the constant mode, and the second one fixes the hourglass mode to be 
el imin ated later. Let us fix u* 1 equal to zero at the origin and at 
the next point on the boundary (coordinates : h,0) (Figure 18. a). 
According to the error estimates (3.30) and (3.31), we may write 

u h «= ^ + XH + 0(h 2-S ) (5.1) 

and therefore, if we normalize H such that its nodal values are 0 or 
1, X measures precisely the value of u h at (h,0) (Figure 18. b), and 
approaches u(h,0) 

X = u(h, 0) + 0(h G ) (5.2) 

2 

But u(h,0) is 0(h ) for a smooth enough solution (u(0,0) = 0 , 
8u/3 n (0,0) = 0) and using L°°-estimates 12 , 0 can be evaluated to 
2-s, e arbitrary. Finally, we have the estimate 

2-e 

X = 0(h ), e arbitrary (5.3) 


Also, the choice of H leads to 



and therefore we obtain 


|| xh|| s = 0(h 2-S " E ) s = 0,1 


(5.5) 


and that proves that the post processor contribution XH can be neg- 
lected if the fixed nodes are chosen as indicated for this type of 
boundary condition* The error estimates of Theorem I still hold up to 
within h G * 

Unfortunately this remark has two major drawbacks: it supposes 

2 

that u is smooth (u G H /(ft)) and it is not valid to 9-node elements 
that will later be discussed. 

Before discussing the implementation of (3.30), we indicate that 
this projection can be simplified. Indeed, taking v* 1 = H in (4.25) 
we obtain 


(f i H >0 

IIhII* 


H II <|[ f h || i^° 


— II ^11 

24 11 U 0 


(5.6) 


(f^ 1 H) , || H |I _ , , 

1 ° H H, < ||f h || 0 — °= H^ll, 


II H || ^ 


II h || t iiT 


Therefore we have 


(f h 

II H|l s < C|| f h |l 0 


2-s 


II H || ^ 


s = 0,1 


(5.8) 


and this term can be neglected without affecting the estimate of Theorem 


VIII. The formula used in the post processor is then 


u* 1 = u* 1 - XH (5.9) 

X = (u* 1 , H) 1 || H I] ~ 2 (5.10) 


In order to preserve the efficiency of the method provided by 
under integration, one must find an efficient way to compute the para- 
meter X in (5.9). One way that suggests itself is to calculate the 
inner products of (5.10) using numerical integration. The use of a 
one point rule would be absurd and would lead to a ratio 0/0. The use 
of a 4 Gauss point rule has been numerically implemented and gives good 
results (similar to those to be presented next) but cost of this inte- 
gration is expensive, as shown in Table 3. This method is therefore 
rej ected. 

We shall now describe a more efficient method with related numer- 
ical results shown in the next subsection. This method relies on the 
fact that, for the bilinear element, the stiffness matrix can be decom- 
posed into two parts, one of which contains H in its kernel. The 
other part is such that the image of H is cheap to calculate. 

This decomposition proved in Section 2.2.2 can be written as 


R exact _ ^under + - T 

e e e -~e ie 


(5.11) 


exact under 

where K and K respectively are the exact element stiffness 

~e ~e 

matrix and its under-integrated form, and are obtained from 


(2.27) and (2.28). In particular 



Si 



s 2 


(5.12) 
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Even though given by (2,27) is still difficult to calculate, note 

that, during the element calculations, the vectors b^ and b^ and the 
Jacobian \& e \ are necessarily computed. Therefore y^ is very easy 
to calculate. 

The inner product (u* 1 , H)^ may, therefore, be calculated by 

exact 

using the decomposition (2.1) of K . Introducing the nodal vec- 

tors U and H associated with the function u* 1 and H , we have 


(u 1 , H) = U T K exact H = U T ‘ Ze Y ‘Y T * H 
1 ~ ~ - e eie ie 

= El(y T • U) (y! • H ) 
e ~e -e ~e -e 


where U and H are the values of U and H at the nodes of the 
~e ~e 

element e . We note that if the values of H are 4-1 or -1 , the 

T 

scalar vector product y^ is always ±4 • Therefore 


(H* 1 , H) = 4Z ± e Z y 1 u 1 
1 e . 'e e 

e 1=1 


(5.13) 


(H, H). - 16Z e 
1 0 e 


(5.14) 


These expressions are still exact since no approximation has been 
made on . If we suppose that the Jacobian of the element is approx- 
imately constant (true for parallelogram element), is simply ex- 

pressed as 


e 12|ftT ( ^L h + &2 ~2> 


(5.15) 


The calculation of the approximate projection can be summarized in 



the following algorithm: 


• Loop on Elements 

A Calculate y , e using (2.2) and (2.6) 
~e e 

T — 

Calculate ±e (y • U ) 
e e e 


J-Add 


{ 


K ± • V J 

1 e e e 


*2 ' X 2 + E e 


\ - X 1 /4X Z 


Loop on Nodes 


t 


u 


‘Node 


-hi 
= u 


Node ± 


X 


Remark : The notations previously used are essentially those found 

in the work of Belytschko and co-workers [4,5] on stabilization methods. 

These methods rely on the decomposition (5.11) but the stabilization 
T 

term ey*y is a-priori added to the under-integrated matrix to prevent 
the spurious modes from the kernel of the stiffness matrix; whereas our 
control method uses the very same term a-posteriori, after solving with 
the under integrated matrix. Therefore, our method seems to be cheaper 
than the stabilization methods as summarized in Table 4. 

2.5.2. Numerical results . 

2. 5. 2. a. Regular mesh of 4-node elements . In order to illus- 
trate what has been stated, we have considered the Laplacian problem 
solved on a square domain partitioned into N^(=h ^) subdomains, for 
various values of N and we have studied the norms of the difference 
between the solution obtained with a full integration u* 1 (4 point rule 



and with under integration (1 point rule) u* 1 and u^ (before and after 
post processing). The results are shown as plot of Log|| u -u || or 
Log|| u-u h || g in function of |Log h| , for s * 0,1 . Data of vari- 
ous regularities have been used: 

i) f^ is a C^-f unction, but not : 

r f^x.y) = f(l-x) - -y y if Y 1 (x,y) > 0 

Lf^x.y) - y(f (i-x) - y if Y i (x ’ y) - 0 

where the (^-discontinuity line is 
(x,y) = y (1-x) - y 


ii) f^ is a non-continuous function 


f 2 (x,y) = 1 

if Y^x.y) > 0 

f 2 (x,y) = -2 

if Y^x.y) < 0 


where is the same as in i) . 

Remark : Both of these functions satisfy the compatibility 

condition (3.2). 

Results obtained with the continuous function f^ are shown in 

Figure 19. When the solution has been treated by the post processor 

2 1 

(Fig. 19a.), bor both L and H norms, the representing points lie 

on lines of slope 2. This proves that whereas the estimate (1-13) is 

2 1 
optimal for the L norm (s = 0), it is not in the H norm and seems 

to be in fact better than what was expected in our study. This does 

not affect in any case the comparison with the exact solution (1.14). 

Figure 19. b shows the comparison with the crude solution u^ 1 , 

obtained with two fixed nodes, and not treated by the post-processor. 



Slopes 2 and 1 are observed and the loss (1. instead of 2. for the 

norm) corroborates the final remark of Section 3.32. 

When the function f^ is used (Figure 20), the points show os- 

2 

cillations around two lines of slope 2. (for the L -norm) and 1.65 
(for the H^-norm) proving that the estimate (3.31) still holds (Fig. 

20. a). When the solution has not been treated by the post-processor 
(Fig. 20. b), the slope 1.65 becomes 1. as expected. 

The next series of examples was intended to study the influence 
of a singularity (at the origin) for a unit square domain regularly 
partitioned. The data functions are of the form 

f a (x,y) = r a - C , a> - 2 (5.16) 

where C is a real number chosen such that the equilibrium condition 
(2.2) is satisfied. The family {f^} represents various regularities 
of data: 

f 6 H s (a) a > s - 1 (5.17) 

The result shown in Fig. 21. a is a plot of a (regularity) versus a 

h h 

(rate of convergence of || u - u || ft -) . The pattern of the (a, a) 

s u f J- 

plot seems to show a linear increase of slope 1 towards the maximum 

2 2 

value 2 reached for f G L (a“ “1) for the L -norm (s=0) . As far as 

the H^-norm (s=l) of the error is concerned, the linear increase of 

2 

slope 1 reached 1 for f G L but keeps increasing towards 2 . This 
shows that the expected error estimate 


||u h -s h !L ich k IU II 

s — m 


9 


s = 0,1 


(5.18) 



where 


k - 1 + min (1, m) - s (5.19) 

is not optimal for $ = 1 and m > 0 . The estimate (3.32) remains 
optimal however. 

In conclusion, these numerical results prove that the method is 
accurate for regular mesh and that no accuracy is lost. 

2.5*2.b. Regular mesh of 8- and 9-node elements . Since the be- 
ginning of Part II we have not discussed the underintegration of the 
stiffness matrix of the 8-node elements. It is well known that this 
matrix is not rank-deficient, and the practice of the under integration 
has been widely used with good results when the mesh is regular. Since 
there is not any spurious mode, the a-posteriori control previously 
described is not needed. 

Unfortunately, the method of proof presented in Section 2.4 cannot 
be used because this element does not possess the nice tensor product 
properties on which the method relies. The only hope for a proof of 
convergence would be to obtain the result as a by-product of a result 
for the 9-node elements. 

As far as this element is concerned (9-node element) , we have 
proved (Section 2.2.3) that the underintegration of this element leads 
to a rank-deficient matrix; in fact, the procedure described in Section 
2.3, for the resolution of the under integrated problem and the projec- 
tion of its solution is completely applicable to a mesh of 9-node ele- 
ments. Thus, Theorem VII is valid and the projection defined in (3.3) 
can be used to eliminate the spurious mode. As far as the existence of 


a convergence theorem is concerned, one can establish generalizations 
of (4*1 6) and (4.17) to 3-node, one-dimensional elements: there exist 

^ , $i » and such that 


(,h - ♦Vb ■ a i (vh - 

(v h , ^' i ) 1 = tjr) 0 


Vv h 


Unfortunately, the basis functions and ^ are different for the 

2 1 

L -under integrated and H -norms and a lemma as Lemma 4.1 cannot be ob- 
tained. 

However, in this subsection we will show numerical results 

obtained by use of the projection (3.30) for regular meshes of 9-node 

elements. Note that two types of control have been tested with similar 

T 

results: the control only involving the term in y.y predicted by 

T 

Belytschko [3] and the complete control calculated with > i = 7, 

8,9. (See Section 2.2.3). The results obtained with either of them are 
similar for this operator (-A) . 

For 8 and 9-node elements, the optimal rates of convergence are 
given by 


II u-u h || s < C h k II f || m s = 0,1 (5.20) 

where 

k = 2 + min(l,m) - s (5.21) 

3— s 1 

and the best rates of convergence 0(h ) are obtained when f 6 H (ft) 

The results obtained with functions presenting a singularity line (such 
as the functions f^ and f ^ previously defined and others) are pre- 



sented in Table 5 (first and second lines). We obtained 1.99 and 1.74 

2 

for a discontinuous function (f G L ) , then 2.43 and 1.97 for a con- 

1 1 12 
tinuous, not C , function (f G H ), 2.95 and 1.95 for a C , not C 

2 00 

function (f G H ) and finally 4 and 3 for a C data. Therefore, the 

2 

rates 3 and 2 are reached when f is at least H or equivalent when 
the solution u is in . In this case, the convergence rate (5.1) 

does not seem to be reached. 

The second series of data involving the singularity at the origin 

(5.16) has been tested and results are shown in Fig. 20. b. The pattern 

of the (a, a) plot shows linear increases of slope 1, the predicted 

values 3 and 2 are reached for f G H^-(r) according to (5.21), but the 

2 

maximum values 4 and 3 are reached for f G H (ft) . 

2.5.2.C. Irregular mesh of 4- and 9- node elements . Finally, the 
method has been tested on the quarter unit disk shown in Fig. 22 with 


f 



a > - 2 


The plot (a, a) is shown in Fig. 23 and we can point out: 

• The general pattern is respected (linear increase towards a 
maximum value) 

• The maximum values 4 and 3 (9-node elements) are reduced to 
values slightly lower than 3 and 2 • 


2.6 Excitation of Spurious Modes 

The previous sections were devoted to the study of the Laplace 
equations with Neumann boundary conditions. The choice of these bound- 
ary conditions is convenient for the analysis of the hourglass instabil- 



ities because these modes appear explicitly in the kernel of the under- 
integrated discrete operator* When Dirichlet conditions are applied 
on a part of the boundary, even though the kernel of the underintegrated 
stiffness matrix is not rank-deficient, instabilities may appear* 

In this section we would like to study the influence the boundary 
conditions have on the solution of the underintegrated problem, and 
obtain results analogous to Theorem VIII* Also we would like to explain 
how the oscillations may be excited in certain problems. The method of 
proof Is similar to that presented in Section 2.4. For various boundary 
conditions, we are able to exhibit the exact eigenvalues and eigen- 
functions of the various linear and bilinear form involved. The expla- 
nation of the excitation of oscillations will result from the comparison 
of these eigenvalues. The procedure also allows us to study the under- 
integration of the operator -A+l and the control of resulting spurious 
modes. Numerical results will illustrate the theory. 

2.6.1. The underintegrated problem with Dirichlet or mixed bound- 
ary conditions . This section is devoted to a generalization of the 
results obtained in Section 2.4 to the Laplacian equation with Dirichlet 
or mixed Dirichlet-Neumann boundary conditions. Only proofs for the 
Dirichlet case will be given in this section, but their equivalent for 
the mixed case can be found in Appendix A. 

The Dirichlet case is simpler than the Neumann case because the 
hourglass mode does not belong to the new approximation space defined 
to handle the boundary condition. Therefore the stiffness matrix is no 
longer singular and can be normally inverted. In the variational for- 
mulation, similar to (3.7), the projection of the data function is not 



necessary and the problem P 


is written: 


where 


u 

(P ) Find u € such that 


. v\ >h = (f h , v h ) 0>h 


Vq = {v h /v h 6 C°(n) , v h | nij 6 Q 1 (fiij) 


v Un = 0} 


V h 6 vjj ( 6 . 1 ) 

, 1 < ±,j < N , 


Remarks 

i) The fact that (*,*)^ ^ Is not singular on vjj does not make 
the problem classically eliptic in the sense that the constant in the 
Lax Milgram Theorem is h-dependent. 

ii) When Dirichlet or mixed boundary conditions are applied. 


Ker « Ker ^ = {0} 

Thus the post processor is not justified anymore and we will compare 
directly u* 1 and u* 1 . 

This comparison is carried the same way as in Section 2.4 and a 

basis of the approximation-space can be obtained. One useful 

1 2 

basis is the common eigen-basis of the matrices of the H -, L -, and 
12 

underintegrated H - or L - inner product. Let us consider the N-l x 
N-l matrix 



D(k) = 


0 


-1 


-1 

2k 



The values for which det D(k) vanishes are 


i i7r 

k i " cos TT 


l<i<N-l 


( 6 . 2 ) 


and the corresponding vectors (D(k^)v^ - 0) are 


Vi^sin^} 


(6.3) 


j=l, N-l 

Let 4>1 = 4> ± (x) , x G £p,lj , denote the piecewise linear function 


associated with the vector : 


where 


4> 1 (jh) © pin ^ 1 < i,j < N-l 

spanU 1 }^^ - vj >0 


V^ Q = (v h G C°(0,1), v h (0) = v h (l) = 0 


(6.4) 

(6.5) 


v is linear on eh, (e+l)h , 0<e<N-l} 


From this point, the remainder of the proof goes as in Section 2.4 and 
the variational problem and its underintegrated formulation can be ex- 
plicitly solved and the decomposition (4.29), (4.30) and (4,31), (4.32) 
are obtained for 1 < i,j <_ N-l , and we finally obtain the result for 
Dirichlet boundary conditions: 

2 

THEOREM IX : Let f be a function in L (Q) . Let u be the 

solution of P : 

P : Find u G / (u,v) ][ = (f,v) Q 

h 2 h 

Let f be the L -projection of f onto V 


Vv G hJ 


( 6 . 6 ) 


and let u be the solu- 



P* 1 : Find u^S vj/ (u*\ v h ) 1>h = (fh » vh) o,h 6 V Q (6,7) 


Then we have the following error estimate: 


— h j 
u-u 


l 8 — 


< C h 2 - s Ilf || 0 


s = 0,1 


□ 


( 6 . 8 ) 


This theorem proves that the use of the underintegrated matrix does not 
affect the rate of convergence of the solution. The method is there- 
fore accurate and efficient. 

Varuous regularities of data have been tested for meshes of 4-, 8-, 
and 9-node elements, with various boundary conditions. Results are 

summarized in Tables 5 and 6. They indicate that the optimal rates of 

2 o 

convergence for f 6 L (£2) for the 4-node case and f 6 H (Q) for the 

8- and 9-node case. 

2.6.2. The under integration of the operator -A+l . In this sub- 
section we consider the underintegration of the operator associated with 
the problem 

Pq : Find u 6 H^(£2) such that 


- Au.+ u = f in ft 


_3u 

3n 


on 


3£2 


(6.9) 


The usual variational formulation of P^ is 


P : Find u 6 H (£2) such that 


(U, v) L + (u, v) 0 = (f, v) 0 , V 6 H 1 ^) (6 . 10) 


The results of existence, uniqueness of P are well-known and so are 
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"the ones for its discrete formulation: 


: Find u* 1 G V* 1 such that 

f h h v , , h h v h v \ / h h 

(u , V ) 1 + (u , V ) 0 = (f, V ) 0 , Vv G V (6.11) 

1 

where V is an approximation of H (ft) using bilinear elements. 

The under integration of (•/> 1 + (•,*) Q leads to the following under- 
integrated problem: 


: Find U* 1 G V* 1 such that 


i — h h — h h N 

(u , v ) 1>h + (u , v ) 0>h 


where the choice of approximation space 

— h h 

\T = V/H 


(f ’ vh) 0,h ’ ^GV 1 


( 6 . 12 ) 


(6.13) 


is justified by 
h 


(v\ H) x h + (v, H) Q h = 0 , V v h 6 V h 


(6.14) 


Then, the method of proof used in Section 2.3 allows us to obtain the 
existence and uniqueness of u 11 . A projection similar to (3.30) can 
be obtained by analogy: we have 


(u h , H) l + (u h , H) Q = (f, H) Q (6.15) 

we therefore construct the projection as: 

u h = mi* 1 * u + X Q H (6.16) 

(u h , H) x + (u\ H) q - (f, H) q 


(6.17) 



This defines uniquely X^ as 


X_ = 


(f,H) 0 - (u'.H^ - (u,H) 0 


Ml + IIhIIJ 


(6.18) 


Similarly to what was done in Section 2. 5.1, we can use (4.25), (5.6) 
through (5.8), simplify X Q without any loss of accuracy and still use 
(5.9), (5.10) for the projection 


~h — h , 
u = u - AH 


X = H) 1 II H 11 “ 


-2 


(5.9) repeat 
(5.10) repeat 


The proof of the convergence of u towards u is again done by 

direct calculation of u* 1 and u* 1 : the explicit resolution of 

and (P. + it) leads to : 
h 


i J 1 + S. + Bj ij 


0 < i,j < N 


(6.19) 


and 


r ~ 
u 


a. a. 

1 3 


i»j 


°i°j + “i B j + a j P i 


f ii 


0 < i,j £ N 
(i,j)*(N,N) 


^ U NN U NN 


( 6 . 20 ) 


h h 

These decompositions allow us to obtain u - u as done in Section 

2 

(2.4). Provided that f € L (ft) we can obtain 


l|u h -G h || s <Ch z S ||f || 0 s = 0,1 


2-s 


( 6 . 21 ) 


Once again, the under integration does not seem to affect the rate 
of convergence. The result can also be obtained with various boundary 
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conditions. Numerical results are joined in Tables 5 and 6 and assert 
the theory. 

2.6.3. Excitation of oscillations . The existence of spurious 
oscillations when underintegration is used is not only encountered 
when Neumann boundary conditions are applied on the whole boundary. In 
this subsection, we would like to analyze precisely how modes that os- 
cillate with wavelength of order h are excited when under integration 
is used, whereas they are damped when the integration is exact. 

For this discussion we consider the unit square 

Q Oo.l [xj0,l[ (6.22) 

discretized into NxN elements. We consider the Laplace equation on S2 

-Au = f in J2 

u = 0 on A {x = 0} 

3u = g on 3fi/{x = 0} 
dn 

For the first time we include two kinds of load: body forces -and sur- 

face loads, and we will observe separately the effects of each of them. 

The eigenfunctions associated with these particular mixed boundary 
conditions are constructed as in Section 2.4. 

X^Cx.y) « ^ i (x)^(y), l<i<N (6.24) 

0<j<K 

where <J>^ is defined in (4.15) (associated with Neumann boundary con- 
ditions at both ends) and is similarly defined (see Appendix A) . 

These functions are defined through sine and cosine functions and there- 
fore oscillate. Among them we will distinguish "smooth" modes with 


(6.23) 
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longer wavelengths (0(1)) from "irregular" modes with shorter wave- 
lengths (0(h)). Smooth (respectively irregular) modes correspond to 
smaller (respectively longer) values of i or j . Examples of each 
extreme are shown in Fig. 25 for N = 10 . 

The resolution of the fully integrated problem leads to the 

search for coefficients u . . such that 

ij 

u h = I u x lj (6.25) 

l<i<N J 
0<j<_N 

The basis {y is an eigenbasis for (*,*)^ and therefore we have 
U ij = A ij [ (f ’ X \ + < ‘ g ’ X J ^ 0 » (6.26) 


where 


with 




i /i7T 7T . 

1 - c ° s( ir - ^ 

o , , ITT 7T x 

2 + C0S( "N" ' 2N> 



1 - COS 0^0 

2 + cos(^) 



(6.27) 


(6.28) 


The values . have been calculated exactly with these formulae 
and their values are reported in Table 7. a for N ■ 10 . The 20 highest 
values are in the shaded zone. We clearly can pbserve that 

i) these values range from the highest value to 1% of this value, 
ii) these values are associated with smooth modes (tensor products 
of smooth modes) . 
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Contrarily, the eigenvalues of irregular modes are smaller and because 
of this, these modes will be damped; only smooth modes will contribute 
in (6.25). 

When the under integration is used, and when g is zero, the solu- 

„ . ~h 
tion u is 


with 


where 


and 



(6.29) 


(6.30) 


(6.31) 


(6.32) 


Again, the values of have been calculated exactly and they 

are reported in Table 7b. The 20 highest values are in the shaded 
zone. The comparison between Tables 7a and b shows that these 20 values 
are approximately the same and they are associated with the same smooth 
modes. In this case, irregular modes will still be damped , and one 
can predict that no oscillation will occur. 

When a load is only applied on the boundary (f = 0, g ^ 0) , 



is now 


u 


±i 


u ij = A ij (s ’ X 


(6.33) 


where 


A ij " + aJ&j 


(6.34) 


Again the values of Ay are reported in Table 7 .c and the 20 highest 
values are in the shaded zone* Among these 20 values, three correspond 
to very irregular modes* In particular, the third value is associated 
with ^10,10 ^ Therefore, we can predict a strong contribution of 
irregular inodes within the solution u* 1 , which will show oscillations* 
Finally, one could wonder if the calculation of the boundary 
integral can be calculated such that (g, X^)q is damped for large 

i and j . Unfortunately, no precise method has been obtained. In par- 
ticular, if the load g is a concentrated load at (x^, y Q ) , then 


(g, X ij ) 0>3Q = X i 3 (x 0 , y 0 ) 

and this value is not necessarily zero. The procedure, consisting of 
splitting g between neighboring nodes, seems to give satisfactory 
results, but is more ad hoc than general. 

2.7. The Practice of Under integration in Linear Elasticity 

We devote this section to the discussion and the effects of the 
under integration of the linear elasticity operator. Our goal is: 1) 
to exhibit the kernel of the under integrated operator; 2) to obtain a 
post-processor formula similar to (3.30) to control, a-posteriori, the 



94 


spurious modes; 3) to indicate how to implement this control, and 4) 
finally, to show numerical results. This study is entirely qualitative 
- the basis function obtained in Section 2.4 cannot be used at this 
point to obtain basis functions for the elasticity operator. However, 
both 4- and 9-node elements will be discussed. 


2.7.1. The kernel of the discrete underintegrated linear elasti- 


city operator. We consider the linear elasticity operator defined by 


a = § T c b 


(7.1) 


where 



(7.2) 


and C is a 3x3 symmetric matrix. In the plane strain case we may 
particularize C : 


( A+2ji A 

A A+2p 
0 0 

In order to exhibit the spurious modes we consider the operator A 

associated with Neumann boundary conditions. In that case, the kernel 

of A consists of the usual 2-dimensional rigid body modes denoted by 

t , t and r 
-x -y 

RBM = span {t^ = (*) ; t * (J) ; r « (“£)} (7.4) 



We consider the problem 
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P : Find u = ( 1 ) G [h 1 ^)] 2 /RBM such that 
~ U 2 ^ ^ 

’ T T f T 

u 8 C 8 v dxdy ■ f • v dxdy 
in ~ Z Z Z ~ IcT 


(7.5) 


where f = (f^, f is a force satisfying the equilibrium conditions 


f G RBM 


(7.6) 


or equivalently 


f dxdy = f. dxdy = 0 

■ q 1 h 

(f y - f x) dxdy = 0 
Jo z 


(7.7) 


The existence and uniqueness of a solution for P are well known. The 
construction of finite element approximations of (7.5) involves the cal- 
culation of the (2Nx2N) stiffness matrix K for a typical element 

~e 

ft , which is given by the formyla 
e 


K =[ N T B T C 3 N dxdy 
* e Jo ~ ~ ~ 


(7.8) 


where N is a vector representing the bilinear (N=4) or biquadratic 
(N=9) shape functions in each element , 1 e <_ E . In computa- 
tional applications is evaluated using an integration rule: 


K = Z W B aT C B a 
= e a=l a 55 * * 


(7.9) 


where, similar to (2*36), 


h‘ 0 -2 

0 ~2 *5 


(7.10) 





and w is the weight at the integration point a . Simple rank con- 

ex 

siderations allow us to predict the rank of . Indeed, since 


rk(A • B) < max(rk A, rk B) 


and rk(A + B) < rk A + rk B 


we have 


rk K < 3 L 
~e — 


(7.11) 


(7.12) 


When the full integration is used, (7.12) does not tell us anything, 
but we know that has the correct kernel containing only rigid 

body modes. However, when underintegration (L=l) is used on 4-node 
elements, we have 


rk T < 3 (7.13) 

Therefore, the 8x8 matrix possesses at least two spurious modes. 

In fact, two is the exact number. Similarly, when under integration 
(L=4) is used and 8- or 9-node elements, we have 

rk K e £ 12 (7.14) 

This inequality predicts one spurious mode for the 16x16 matrix associ- 
ated with 8-node elements, but when the procedure is repeated for two 
neighboring 8-node elements, the spurious modes can no longer exist in 
the global matrix Q»6] . We can also interpret this elimination of the 
spurious mode by noticing that neighboring element cannot share the 
mode ^13^ • 

As far as 9-node elements are concerned, the inequality (7.14) 
tells us that the 18x18 stiffness matrix has at least three spurious 
modes. In fact, there are exactly three such modes and they can be 



shared by adjacent elements. Next the modes will be explicitly de- 
scribed. 


H 


2. 7.1. a. The spurious modes for 4-node elements , 
be the two hourglass vectors defined as 


Let H and 
~x 


H = 
-x 




(7.15) 


* 

where h is the hourglass nodal displacement defined in (2.14). Then 
when L=1 , we obtain from (2.18) and (2.20) 



and similarly 


Therefore, 


B • H = 0 

* ~y 


K' 1 ^ • H « • H = 0 

~e —x ~e ~y 


(7.16) 


These element displacements can be put together to obtain' two 

global spurious modes, also denoted by fl and H and we have : 

~x ~y 


Ker K un ^ er = {span t,t f r, H,H} 

~x -y ~ ~x -y 


(7.17) 


This defines entirely the kernel of the underintegrated matrix and the 
spurious modes for 4-node elements. 

Remark : In problems where symmetry is used for simplifications, 

the kernel of jr un£ * er muS £ respect the symmetry. If one axis of sym- 


* In this section, nodal values and associated functions will be denoted 
by the same letter, the underlining 'V differentiating them. The 
nodal values are expressed component by component. 


metry (say, the x-axis) exists, then 


Ker £ unc * er = S p an (t , H } 
» ~x -x 


If the problem has two axes of symmetry (x- and y-axis) 

Ker K undcr - {0} 

The spurious modes are eliminated by the symmetry conditions. 

2.7.1.b. The spurious inodes for 9-node elements. Let H and 

— ~x 

Hy be the two vectors defined as 


?x = <£ 




(7.18) 


where h is the spurious mode of 9-node elements defined in (2.37). 
Using (2.36), (2.39) and (7.9), we easily get 


(4) (4) 

K .H = IC (H = 0 
~e ~x -e ~y 


(7.19) 


Therefore, H and H are two out of the three spurious modes of 
~x ~y 

(4) 

K g . We remark that the pattern (2.37) defining them does not depend 
on the geometry of the mesh. As far as the third spurious mode denoted 
by 


w. 


W = („ x ) 


w„ 


(7.20) 


is concerned, one can show that the equations defining it are 


K aT _ . aT 

Si * ~i " h ' Z2 zi ' Z2 ■ Z2 


= b? T . w^ + b“ T « w - 0 ; a=l,4 (7.21) 


or equivalently : 



Note that for this system of 12 equations, we have 18 unknowns. If 

we add 5 orthogonality equations between W and t , t , r, H and H , 

— —x ~y ~ —x -y 

the system will define only one W (up to within a multiplicative fac- 
tor) . For a general geometry of , one cannot exhibit an explicit 

form for W ; however, when is a quadrilateral, and when x and 

y are of the form 

x = (x^ x 2> x 3> x^, ^(x 1 + x 2 ), %(x 2 + x 3 ), 

J$(x 3 + x^,), %(x^ + x^) , + x 2 + x 3 + x^)) (7.23) 


we can prove that one candidate for W can be written as 





w 2 = -T x’ J 


(7.24) 



and x T , y T are the vectors constructed with the first four components 
of x and £ . An example on W for a geometry satisfying (7.23) is 
shown in Figure 25 and can be constructed as follows: 

i) the displacement of a mid-side node is normal to the side. 


alternatively inwardly and outwardly oriented, with magnitude 
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proportional to the length of the side. 

ii) the displacement of a corner is obtained by multiplication 
of -2 of the sum of the two displacements of the closest 
mid-side nodes. 

iii) the displacement of the centroid is zero. 

On a square, the pattern of W is well-known: 

/- 2 , 2 , 2 , - 2 , 0 , - 1 , 0 , 1 , 0 \ 

W- (7.26) 

\ 2 , 2 , - 2 , - 2 , - 1 , 0 , 1 , 0 , 0 / 

Contrary to 8-node elements, and because of the presence of H 
and H y , this mode can "propagate" from one element to another. For 
example, on a square mesh, if the nodal displacement vector is W 
given by (7 .26) on an element fig , then the displacement vectors W 
+ 3H + t and W - 3H y - t on the elements to the right of fi 

and above JJq allow us to construct a continuous global displacement 
also denoted W , on the mesh as shown in Figure 26. 

We finally have 

Ker k 1 ™* 161 = span t , r, H £ , H y , W} (7.27) 

Remark : Similar to what we have with 4-node elements, the exis- 

tence of one axis of symmetry (say, the x-axis) reduces the kernel of 
the underintegrated stiffness matrix: 

under 

Ker K = span (t^, H 2 + H 3 ) (7.28) 

where 

i - 3/2 ( 5x - 
H 2 - -3/2 (H y - t y ) 




101 


~ H 3 = ? + 2(t x - V ' 


(7.29) 


have been chosen such that the displacements of these modes are zero at 
the intersection of both axes for a square mesh* Contrary to the 4-node 
case, we still have a spurious mode when two axes of symmetry exist: 


Ker K un< * er = span {H^ + + H^) 


(7.30) 


This mode is shown in Figure 27. 

It is also important to point out that whereas the pattern of 

the spurious modes and are independent of both the geometry 

and the element, the mode W depends upon both of them. Moreover, we 

can see by construction on a square mesh that the amplitude varies 

strongly when we consider successive elements. In fact, the pattern we 

may observe is a succession of pattern H and H with increasing 
. — x ~y 

amplitude. 

2.7.2. The a-posteriori control in linear elasticity . In this 
subsection we wish to generalize (3.30) with regard to the discrete 
operators, using various kernels discussed in the previous subsection. 
We consider the general case where 

Ker K UtlC * er = rbM © span{H^, i * l,l} (7.31) 


where I may have the values 1, 2 or 3 . We recall that for 


I = 1 , we obtained a control formula similar to 


_h -h v 

u " u " aO^, H 1 ) 


H, 


(7.32) 


where the bilinear form a(*,») was obtained in the variational formu- 
lation of the initial problem. This projection satisfies: 
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a(G h , H x ) = 0 (7.33) 

or, in other words, u* 1 is orthogonal to the spurious mode. We gener- 
alize this property to the elasticity problem by supposing the projection 
to be orthogonal to all the spurious modes. Therefore the control will 
consist of looking for I constants X^ (i - 1,1) such that 


r Q h = ^ - i x . h 
J - - 1-1,1 

l a(u h , = 0 for i-1,3 


(7.34) 


This leads to the system of I equations with I unknowns : 


Find X i 


Z 

J-l.I 


, i * 1,1 such that 
Xj a(H ± , H ) = a(u^ , 


H i> 


i=l,I 


(7.35) 


The computations involved in the control are computations of products 
of u* 1 and the spurious modes by themselves. The implementation of 
these computations are to be discussed in the next section. 

2.7.3. Implementation of the spurious modes control . For the 
computation of the coefficients in (7.35) we again use the decomposition 


K f “H - K U ” der + K> 
where K under satisfies 


then 


under 

K 



aCu* 1 , H^) 


— T * 

l U * K • H 
e=l,E 


(7.36) 


(7.37) 


(7.39) 



The expressions used for are next given for 4- or 9-node 

elements, 

2, 7. 3. a. Control for 4-node elements . For the operator defined 


in (7.1) and (7.2) 

with 


N 

C 12 

C 13 

2“ C 21 

C 22 

C 23 

\ C 3! 

C 32 

C 33 


(7,40) 


we have the exact decomposition for any geometry of £2 


exact _ under 

& — ■ IX + 

-e - 


a n VI 


a 21 VI 


a i2 VI 


a 22 VI 


where 


Hi 


C ll*«^ C 13^31 ) V^V^'ix +<C 12 4C 33 )e x y 4<: 32V 

C 31 £ xx +(C 21 +C 33 )e xy' K: 23 e yy ;C 33 e xx +(C 32 4C 23 )e xy +C 22 e yy 


(7.42) 

The vector y and the F's are defined in Section 2 ((2.41) and 
(2.50)). For practival use, the expressions (2.51) are used for e . 
For linear isotropic linear material, C is given by (7.3) and 




/(X+2u)e + y e 

1 xx yy 


v e. 


m e 


xy 


xy 


P e +(A+2P)e 

xx ' yy/ 


(7.43) 


This expression of 
relationship: 


H 


can be compared to the general strain-stress 
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M 


(X+2y)e + ye 

x y 


y e 


V e 


xy 


xy 


ye +(X+2y)e 
x y 


(7.44) 


An algorithm similar to the one presented in Section 2.5.1 can be con- 
structed. It involves the computation of y , £ and a , then the 
computation of a(H^, H^) and aCu* 1 , H^) , and finally the coeffi- 
cients X^ are obtained by resolution of a NxN system, N measuring 
the rank deficiency of ^ under (n= 1 or 2). 

2.7.4,b. Control for 9-node elements . In this subsection, devoted 
to 9-node elements, we would like first to show why the results obtained 
by Belytschko are not sufficient to obtain a generalization of the 
Linear Elasticity Problem, and then to propose an implementation of the 
control that leads to a stable solution converging to the exact solu- 
tion with the optimal rate of convergence. However, for 9-node elements, 
we have not yet been able to obtain a computationally easy way to exhi- 
bit the third spurious mode, and the proposed results are only applica- 
ble to regular discretizations of a domain. 

As far as the stabilization method proposed in ^3 _} is concerned, 
algebra similar to that in Subsection 2. 7. 3. a leads to (7.41) where y 

was defined in 2.41. But, whereas the stabilization matrix constructed 

T 

with the submatrix y*y eliminates H and H from the kernel of 

l l ~x ~y 

the stiffness matrix, it does not take W into account. Indeed, we 
have 


a ii VI 


a 2i VI 


V w - 0 


a i2 VI 


a n v~ 



Therefore this procedure cannot be used to control W . 

In order to obtain an accurate control, we have to consider a 
generalization of (2,43). Now we have 


a s *s 
11 Z± Z± 


a !2 


a 01 s. s. 

21 -l ~i 


a 22 !i f i ' 


• H. t 0 
~3 


for 7 l i l 9 

!<j<3 


where the vectors are defined in Section 2.2.3. Finally, using 

(2.43) and (2.52), we have 


1«>. k (4) + ^ 

»e -e 135 


( C H+ C 33)SgS 9 ^ C i3 +C 32^ ~9 


^ ~9~9 "~l 


L( C 3i+C 2 3^f9?9 ( - C 22 +C 33 J t9Zs 


C 13tlt7 1 . f C 33~8~8 C 32f8f8 T 


i/Hll ~7 ~7 C 13! 7 !7 1| 

r t? a ^ C o c? 


C 31?7f7‘ Z llil s S~‘ ^ C 22 S W C 22 S W 


\ 

} 


(7.45) 


Similarly, for the 4-node case, the algorithm for the computations 

of the coefficients in (7.35) has been obtained and implemented. Numer- 

T 

ical results agree with our presumptions concerning a "y*Y "-type of 
control and incline in favor of the decomposition (7.45). On a square 
domain discretized with NxN elements, we have calculated and compared 
the solutions obtained with full (exact) and underintegration for vari- 
ous boundary and symmetry conditions. The rates of convergence were 

2 

calculated by comparing the error norms (s=0: L /RBM norm; s=l: energy 

2-s 

norm) obtained with N=5,6 and 7. We consistently got the rate 0(h ) 

T 3 — s 

using a y*y decomposition and 0(h ) with (7.45) for homogeneous 



materials under the action of gravity (f G C ) ; the order 3 -s being 
optimal, we may conclude that the method presented below is accurate. 

It is also efficient: for one second taken for the fully integrated 

stiffness matrix, only .61 are taken when the underintegration is used 
and only .05 seconds are taken for the control. 

Remark : The analysis of the excitation of the spurious modes 

carried for the simple Laplace equation cannot be done for the elasti- 
city because we are not able to exhibit eigenbasis of the discrete 
operator for neither 4- nor 9-node elements. Numerical computations 
a seem to indicate that the same phenomenon occurs: several "irreg- 

ular" modes appear within the smooth, high wavelength modes, and are 
therefore excited. The shape of these modes and their mathematical 
knowledge would allow their elimination or damping. 

2.8 Conclusions and Further Research 

The underintegration seems to be a very attractive way to obtain 
more efficient computations in solid or fluid mechanics. The spurious 
modes this practice introduces and the precise way they are excited have 
long remained unstudied. Several authors previously mentioned proposed 
several interpretations based on the intuition. We have here tried to 
study this phenomenon from a rigorous mathematical point of view, and 
we have precisely answered all the questions concerning one simple prob- 
lem. Unfortunately, the algebra involved in more sophisticated problems 
(9-node elements, linear elasticity) does not allow us such a complete 
study. We would like to indicate that several generalizations of our 
results will help in the control of the spurious modes: i) a discrete 



eigenbasis of the linear elasticity would allow an interpretation of 
the excitation of spurious modes, and hence a way of damping them; 
ii) an accurate control formula for 9-node elements would help in 
a-priori as well as a-posteriori control of the widely known spurious 
mode W • This could also help in preventing bad behavior of 8-node 
elements in certain geometries. 



APPENDIX A 


As far as mixed boundary conditions are concerned, we suppose 
that a Dirichlet boundary condition is applied at 0 and a Neumann 
boundary condition at 1 . For the interval [j),l] , we consider the 
NxN matrix 



The values for which det K(k) vanishes are: 

k i = cos (ff + lf)’ 1 l 1 £® (A. 2) 

and the corresponding vectors (D(ki)vi = 0) are : 

v i = (sin } 1 £ j £ N (A. 3) 

The corresponding approximation space V k 1 with basis {<J>^} is con- 

structed as in £25]or in Section 2. 4* Then, depending upon the sides 

where the various boundary conditions (D or N) are applied, tensor 

product of v! 1 , v! 1 or v!* t are to be considered. The results of 
i 1 , u I 

Theorem II hold for the Mixed Problem. 
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Method 1: 
Method 2: 
Method 3: 


Table 2: Norm- evaluation obtained by: 

Q^/P^ elements 
composite elements 

18/Pj^ elements and filtering of the pressures 
by using only the centroidal value 


Exact Solution: 


IpIIlW* =10 ° 


|-|| = 175.5942 


h = 0.0625 


H p e I' L 2 (fi) /Pv 

H p ~ p e Hl 2 (S2)/TR 

H p_p e 11l 2 (J2)/TR 

l p L 2 cn)/m 

Method 1 

167.1254 

20.0310 

0.1141 

Method 2 

171.5448 

36.3181 

0.2068 

Method 3 

171.5845 

26.6219 

0.1516 
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Table 4: Operations Cost per Element for Both 

Stabilization and A-posteriori Methods 


Stabilization Method 

Operations Cost 

A-posteriori Method 

Computations of Y e 

2 Ox ; 21+ 

2 Ox ; 21+ 

Computations of ^ e *-Y e 

Multiply y e • y e T 

16x 

Ax 

T 

Multiply • Y 0 

Multiply £ e ’Y e Y e T 

16x 

lx 

T 

Multiply e • Y e 

Add K e + c e Y fi Y l 

16+ 

2+ 

Y ± e U T y 
A dd 1 6 6 6 

Y 2 + E e ' 


j 

4+ 

Then: (4 nodes/ element) 
Add ^ ± X 

TOTAL 

52x ; 37+ 

29x ; 27+ 

TOTAL 
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Table 5: Rate of Convergence ] Log |j e ] | |v.s|Log h| (0-0,1) 


for 8- and 9-node Elements 


| REGULARITY 

I 

I BOUNDARY CONDITIONS 
I AND ELEMENT 


f 6 L 2 I f e H 1 | f G H 2 
t H 1 j 0 H 2 ] 6 E 3 


fee 


INEUMANN, 9-NODE 
l+SPECIAL PROJECTION 


1.99 12.43 12.95 

1.741 1.971 1.94 


INEUMANN, 8-NODE EL. 


1.99 12.00 12.97 

1.791 1.931 1.95 


IDIRICHLET, 9-N0DE EL. 


I DIRICHLET, 8-NODE EL. 


[MIXED, 9 -NODE EL. 


I MIXED, 8-NODE EL. 


2.35 12.85 12.99 

1.471 1.991 1.99 

2.30 12.71 12.99 

1.461 2.121 1.99 

2.00 12.00 13.83 

1.671 2.281 2.84 

2.00 12.00 13.84 

1.741 2.271 2.84 
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Table 6: Rate of Convergence [Log || e || |v.s.|Log h| for 

s 

4-node Elements (s=0,l) 


* . * * * * 

I REGULARITY | I ' I „ I 

I IfeL lfGH lfeCl 

I OPERATOR AND I u „1 I . 2 I I 

I BOUNDARY CONDITIONS I * I C I I 

* * * * * * 

I I NEUMANN 11.99 12.00 |2.00 I 

I I +PRO JECTION I 1.611 2.001 2.001 

| - A I DIRICHLET 11.99 12.00 12.00 I 

| | | 1.501 1.851 2.001 

* * * * * * 

I I MIXED 12.00 12.00 12.00 ! 

| | | 1.501 1.991 1.991 

* * * * * * 

| I NEUMANN 12.00 |2.00 12.00 I 

I I +PROJECTION I 1.501 2.00! . 2.001 

* * * * * * 

| - A + 1 I DIRICHLET {2.00 1 2.00 12.00 I 

| | | 1.501 1.851 2.001 

* * * *- -* * 

| | MIXED 12.00 12.00 12.00 I 

| | | 1.501 1.851 2.001 

it * * * * * 
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i 

4 

1 

Table 7 : Arrays of Eigen Values A 

Table 7 a . Array A ^ 

2 3 4 5 6 7 8 

A 

9 

. . and 
ij 

10 

J 

0 1 

40.45 

4 . 42 M 54 0:75 

0.43 

0.27 

0.18 

0.13 

0.10 

0.08 

1 1 

8705 

3.07 

1,34 0.70 

0.41 

0.26 

0.17 

0.12 

0.10 

0.08 

2 1 

2.31 

1.58 

0.95 0 . 57 . 

0.36 

0.24 

0.17 

0.12 

0.09 

0.08 

3 1 

1.02 

0.85 

0.62 0.44 

0.30 

0.21 

0.15 

0.11 

0.09 

0.08 

4 1 

0 ' f 55 

0.49 

0.411 0.32 

0.24 

0.18 

0.13 

0.10 

0.08 

0.07 

5 1 

0.33 

0.31 

0.27 0.23 

0.19 

0.15 

0.12 

0.09 

0.08 

0.07 

6 I 

0.21 

0.21 

0.19 0.17 

0.14 

0.12 

0.10 

0.08 

0.07 

0.06 

7 1 

0.15 

0.14 

0.14 0.12 

0.11 

0.10 

0.08 

0.07 

0.06 

0.05 

8 1 

0.11 

0.11 

0.10 0.10 

0.09 

0.08 

0.07 

0.06 

0.05 

0.05 

9 1 

0.09 

0.09 

0.08 0.08 

0.07 

0.07 

0.06 

0.05 

0.05 

0.04 

10 I 

0.08 

0.08 

0.08 0.07 

0.07 

0.06 

0.06 

0.05 

0.04 

0.04 




Table 7 b . Array 




i 

1 

2 

3 4 

5 

6 

7 

8 

9 

10 


J 

0 

1 



mrm i 

0.18 

0.09 0.04 

0.01 

0.00 

1 

1 

Jjp99 

3.02 

T :27 

0.62 

0.33 

0.18 

0.09 0.04 

0.01 

0.00 

2 

1 

W,2U 

1.53 

0.90 

0.52 

0.30 

0.17 

0.09 0.04 

0.01 

0.00 

3 

1 

0.94 0.79 

0.58 

0.39 

0.25 

0.15 

0.09 0.04 

0.01 

0.00 

4 

1 

" 0.47 

0.43 

0 . 3610.28 

0.20 

0.13 

0.08 0.04 

0.01 

0.00 

5 

1 

0.25 

0.24 

0.21 

0.18 

0.14 
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Figure 7. a) The Qi/Po “ element and b) a mesh path needed to construct the 
auxiliary function <j> with values q e for e such that I + J 
is even, so as to obtain (7.1) and (7.2). 
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Figure 9. ; Mesh with composite-elements. 
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Figure 17 . Geometrical interpretation of the projection u n = iru. 




LOG (ERROR) LOG (ERROR) 




LOG(ERROR) 




• L2— NORM 


143 




e 21 . (a»o) Plot for a square domain 
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